Draft version March 2, 2013 

Preprint typeset using I^-T^X style emulatcapj v. 5/2/11 



CLASH: MASS DISTRIBUTION IN AND AROUND MACS J1206. 2-0847 FROM A FULL CLUSTER LENSING 

ANALYSIS* 

Keiichi Umetsu 1 , Elinor Medezinski 2 , Mario Nonino 3 , Julian Merten 4 , Adi Zitrin 5 , Alberto Molino 6 , Claudio 
Grillo 7 , Mauricio Carrasco 8,9 , Megan Donahue 10 , Andisheh Mahdavi 11 , Dan Coe 12 , Marc Postman 12 , Anton 

KOEKEMOER 12 , NlCOLE CzAKON 13 , JACK SAYERS 13 , TONY MrOCZKOWSKI 14 ' 4 ' 13 , SUNIL GOLWALA 13 , PATRICK M. KOCH 1 , 

Kai-Yang Lin 1 , Sandor M. Molnar 15 , Piero Rosati 9 , Italo Balestra 3 , Amata Mercurio 16 , Marco Scodeggio 17 , 
Andrea Biviano 3 , Timo Anguita 818 , Leopoldo Infante 8 , Gregor Seidel 18 , Irene Sendra 19 , Stephanie Jouvel 20-21 , 
Ole Host 20 , Doron Lemze 2 , Tom Broadhurst 19,22 , Massimo Meneghetti 23 , Leonidas Moustakas 4 , Matthias 

B ARTELMANN r> , NARCISO BeNITEZ 6 , RYCHARD BOUWENS 24 , LARRY BRADLEY 12 , HOLLAND FORD 2 , YOLANDA JlMENEZ-TEJA 6 , 

Daniel Kelson 25 , Ofer Lahav~ , Peter Melchior 26 , John Moustakas 27 , Sara Ogaz 12 , Stella Seitz 28 , and Wei 
O . Zheng 2 

Draft version March 2, 2013 

£3 ■ ABSTRACT 

We derive an accurate mass distribution of the galaxy cluster MACS J1206. 2-0847 (z = 0.439) from 
\Q , a combined weak-lensing distortion, magnification, and strong-lensing analysis of wide-field Subaru 

BVR c I c z' imaging and our recent 16-band Hubble Space Telescope observations taken as part of the 
Cluster Lensing And Supernova survey with Hubble (CLASH) program. We find good agreement 
Q"^ , in the regions of overlap between several weak and strong lensing mass reconstructions using a wide 

■ variety of modeling methods, ensuring consistency. The Subaru data reveal the presence of a surround- 
ing large scale structure with the major axis running approximately north-west south-east (NW-SE), 

, aligned with the cluster and its brightest galaxy shapes, showing elongation with a ~ 2 : 1 axis ratio in 

the plane of the sky. Our full-lensing mass profile exhibits a shallow profile slope d In Y,/d In R ~ — 1 at 
cluster outskirts (R ^ 1 Mpc h" 1 ), whereas the mass distribution excluding the NW-SE excess regions 
S-j ■ steepens further out, well described by the Navarro-Frenk- White form. Assuming a spherical halo, we 

c/3 ■ obtain a virial mass M v i r = (1.1±0.2±0.1) x 10 15 M Q and a halo concentration c v i r = 6.9±1.0±1.2 

(c v ir ~ 5.7 when the central 50kpc/i _1 is excluded), which falls in the range 4 ,$ (c) ;$ 7 of average 
c(M, z) predictions for relaxed clusters from recent A cold dark matter simulations. Our full lensing 
' results are found to be in agreement with X-ray mass measurements where the data overlap, and when 

, combined with Chandra gas mass measurements, yield a cumulative gas mass fraction of 13.71^% 

at 0.7Mpc/i~ 4 (« L7r 2 5oo), a typical value observed for high mass clusters. 

■ Subject headings: cosmology: observations - - dark matter - - galaxies: clusters: individual 
| (MACS J1206. 2-0847) — gravitational lensing: weak — gravitational lensing: 

strong 

"xl" . 
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1. INTRODUCTION 

Clusters of galaxies are the largest self-gravitating 
systems in the universe. These massive clusters con- 
tain rich astrophysical and cosmological information 
about the initial conditions for cosmic structure for- 
mation and assembly of structure over cosmic time. 
Statistical and detailed individual properties of clus- 
ters can therefore provide fundamental constraints on 
models of cosmic structure formation (e.g., Allen et al. 
2004; Vikhlinin et al. 2009), the unknown nature of 
dark matter (DM, hereafter; Markevitch et al. 2004; 
Clowe et al. 2006) , and possible modifications of the law 
of gravity (Narikawa & Yamamoto 2012), complement- 
ing cosmic microwave background, galaxy clustering, and 
Type la supernova observations (Komatsu et al. 2011; 
Percival et al. 2010; Riess et al. 1998). 

Observations of clusters have provided independent 
pieces of empirical evidence for the existence of DM (e.g., 
Zwicky 1959; Markevitch et al. 2004; Clowe et al. 2006; 
Okabe & Umetsu 2008a; Mahdavi et al. 2007a). A prime 
example of this comes from combined X-ray and lensing 
observations of the "Bullet system" , which is understood 
to be the result of a high-speed collision of two cluster 
components occurring approximately in the plane of the 
sky, displaying a prominent bow shock proceeding the 
cool, bullet-like gas subcluster, lying between the two 
distinct clusters (Markevitch et al. 2004). For this sys- 
tem, the bulk of mass is shown to be associated with 
the bimodal distribution of cluster member galaxies, sup- 
porting that DM is effectively collisionlcss as galaxies on 
sub-Mpc scales (Clowe et al. 2006). Such displacements 
between the gas and mass distributions are quite com- 
mon in merging systems, and exhibit a complex vari- 
ety of merging configurations (Okabc & Umetsu 2008b; 
Mahdavi et al. 2007a; Merten et al. 2011a). 

Substantial progress has been made through numeri- 
cal simulation in understanding the formation and struc- 
ture of collisionlcss DM halos in quasi equilibrium, gov- 
erned by nonlinear gravitational growth of cosmic den- 
sity perturbations. In the standard A cold dark mat- 
ter (ACDM) paradigm of hierarchical structure forma- 
tion, cluster-sized DM halos form through successive 
mergers of smaller halos as well as through smooth ac- 
cretion of matter along surrounding filamentary struc- 
tures (Colberg et al. 2000; Shaw et al. 2006; Gao et al. 
2012). In this context, the hierarchical build up of clus- 
ters proceeds in a highly anisotropic configuration where 
infall and merging of matter tend to occur along pref- 
erential directions (Colberg et al. 2005), leading to the 
emergence of the filamentary network of matter, as ob- 
served in large galaxy redshift surveys (e.g., Colless et al. 
2001; Tegmark et al. 2004; Geller et al. 2011). Cluster 
halos are located at dense nodes where the filaments in- 
tersect, generally triaxial reflecting the collisionless na- 
ture of DM, and elongated in the preferential infall di- 
rection of subhalos, namely along surrounding filaments 
(Shaw et al. 2006). 

The internal structure of DM halos constitutes one of 
the most distinct predictions for the CDM paradigm. 
TV-body simulations of collisionless CDM established a 
nearly self-similar form for the spherically-averaged den- 
sity profile (p(r)) of DM halos (Navarro et al. 1997, here- 
after, Navarro- Frenk- White (NFW)) over a wide range 



of halo masses, with some intrinsic variance associated 
with the mass assembly histories of individual halos 
(Jing & Suto 2000; Tasitsiomi et al. 2004; Graham et al. 
2006; Navarro et al. 2010; Gao et al. 2012). The loga- 
rithmic gradient 73dM = —dhip/dlnr of the NFW 
form flattens progressively toward the center, with an 
inner slope flatter than a purely isothermal structure 
(73D = 2) interior to the inner characteristic radius r s 
providing a distinctive, fundamental prediction for the 
empirical form of CDM halos. A useful index of the de- 
gree of concentration is c v i r = r v - lr /r s , which compares 
the virial radius r v j r to r s . Halo concentration is pre- 
dicted to correlate with halo mass since DM halos that 
are more massive collapse later when the mean back- 
ground density of the universe is correspondingly lower 
(Bullock et al. 2001; Zhao et al. 2003; Neto et al. 2007). 
This prediction for the halo c v i r ~M v i r relation and its evo- 
lution has been examined by several independent large 
scale simulations (e.g., Navarro et al. 1997; Bullock et al. 
2001; Neto et al. 2007; Duffy et al. 2008; Klypin et al. 
2011; Bhattacharya et al. 2011), with sufficient detail to 
establish the inherent scatter of this relation around the 
mean, arising from variations in the formation epoch 
of individual halos of given mass (Wechsler et al. 2002; 
Neto et al. 2007; Zhao et al. 2009). 

Galaxy clusters act as powerful gravitational lenses 
(e.g., Bartelmann & Schneider 2001; Umetsu 2010; 
Kneib & Natarajan 2011), providing a direct probe for 
testing these well-defined predictions of halo struc- 
ture because they are expected to have a relatively 
shallow mass profile with a pronounced radial curva- 
ture. A detailed examination of the ACDM predictions 
by cluster lensing has been the focus of our preced- 
ing work (Broadhurst et al. 2005a; Umetsu et al. 2007; 
Broadhurst et al. 2008; Umetsu & Broadhurst 2008; 
Umetsu et al. 2009, 2010, 2011b,a). 

Recent detailed lensing analyses have shown that 
the projected cluster mass profiles constructed from 
combined weak and strong lensing data have a grad- 
ually steepening logarithmic gradient, in agreement 
with the predicted form for the family of collision- 
less CDM halos in virial equilibrium (Gavazzi et al. 
2003; Broadhurst et al. 2005a; Limousin et al. 2007; 
Broadhurst et al. 2008; Umetsu & Broadhurst 2008; 
Newman et al. 2009; Umetsu et al. 2010, 2011b,a; 
Zitrin et al. 2010, 2011c; Oguri et al. 2012; Coe et al. 
2012). Intriguingly, however, some of these results 
reveal a relatively high degree of mass concentration 
in high- mass lensing clusters (e.g., Gavazzi et al. 2003; 
Kneib et al. 2003; Broadhurst et al. 2008; Oguri et al. 
2009; Zitrin et al. 2011c), lying well above the c v i r - 
M v i r relation for cluster-sized halos (c v i r ~4-5 for CDM 
halos with M v - ir /t 10 15 M© in the local universe) pre- 
dicted by the ACDM model, despite careful attempts 
to correct for sizable (~ 50-100%) projection and 
selection biases inherent to lensing by triaxial halos 
(Hennawi et al. 2007; Meneghetti et al. 2010b, 2011). 
The effects of baryons on the total mass profile are gen- 
erally found to only modify cluster concentrations at 
the ~ 10% level (Mead et al. 2010; Duffy et al. 2010), 
although some studies suggest that low mass systems 
(M vir < 5 x 1O 14 M ) may be significantly affected by 
the effects of baryonic cooling (Fedeli 2011; Oguri et al. 



Umetsu et al. 



3 



2012). This apparent overconcentration of lensing 
clusters is also indicated by the generally large Ein- 
stein radii determined from strong-lensing observations 
(Broadhurst & Barkana 2008; Meneghetti et al. 2010a; 
Zitrin et al. 2011a). These lensing results could suggest 
either substantial additional mass projected along the 
line of sight, due partly to halo triaxiality (Oguri et al. 
2005), or an intrinsically higher-than-predicted concen- 
tration of mass; the latter could imply that clusters 
formed earlier than predicted by iV-body simulations of 
the current concordance ACDM cosmology. 

The Cluster Lensing And Supernova survey with Hub- 
ble (CLASH, Postman et al. 2012) 31 has been in progress 
to obtain accurate cluster mass profiles for a sizable sam- 
ple of representative clusters by combining high-quality 
strong- and weak-lensing measurements, in combination 
with the complementary Subaru wide- field imaging (e.g., 
Umetsu et al. 201 lb, a). CLASH is a 524-orbit multi- 
cycle treasury Hubble Space Telescope (HST) program 
to observe 25 clusters of galaxies at 0.18 < z < 0.89, 
each in 16 filters with the Wide Field Camera 3 (WFC3; 
Kimble et al. 2008) and the Advanced Camera for Sur- 
veys (ACS; Ford et al. 2003), ranging from the UV, 
through the optical, to the IR. Importantly, 20 CLASH 
clusters were X-ray selected to be massive and relatively 
relaxed. This selection avoids the strong bias toward 
high concentrations in previously well-studied clusters 
selected for their strong lensing strength, allowing us to 
meaningfully examine the c-M relation over a sufficiently 
wide mass and redshift range for a cluster sample that is 
largely free of lensing bias (Postman et al. 2012). 

In this paper we present a comprehensive weak and 
strong lensing analysis of the X-ray selected CLASH 
cluster MACS J1206.2-0847 (MACS1206, hereafter; see 
Table 1) at z = 0.439 based on the Subaru wide-field 
BVR c I c z' imaging, combined with our recent CLASH 
HST imaging and VLT/VIMOS spectroscopic obser- 
vations presented in Zitrin et al. (2012), who carried 
out a detailed strong-lensing analysis of the cluster. 
MACS1206 is an X-ray luminous cluster (Ebcling et al. 
2009), originally discovered in the Massive Cluster Sur- 
vey (MACS, Ebeling et al. 2001, 2009). Therefore, it 
is an interesting target for detailed lensing analyses to 
compare with well-studied, lensing-selected clusters (e.g., 
Umetsu et al. 2011b,a; Oguri et al. 2009, 2012). 

The paper is organized as follows. In Section 2 we 
briefly summarize the basic theory of cluster weak grav- 
itational lensing. In Section 3, we describe details of the 
full weak-lensing analysis of Subaru observations. In Sec- 
tion 4, we present results from several scmi-indcpcndcnt 
strong-lensing analyses to test the consistency of our 
strong-lens modeling. In Section 5 we derive cluster 
weak-lensing profiles from Subaru data. In Section 6 
we combine our weak-lensing measurements with inner 
strong-lensing based information from CLASH HST ob- 
servations to make a full determination of the cluster 
mass profiles; then, we examine the radial dependence 
of the cluster mass distribution based on the full lens- 
ing analysis. In Section 7, we assess carefully various 
sources of potential systematic uncertainties in the clus- 
ter mass and concentration measurements, and discuss 
our results along with our complementary X-ray and 

31 http:/ /www. stsci.edu/~postman/CLASH 



Sunyaev-Zel'dovich effect (SZE) observations. Finally, 
a summary is given in Section 8. 

Throughout this paper, we use the AB magnitude sys- 
tem, and adopt a concordance ACDM cosmology with 
O ro = 0.3, Q A = 0.7, and H = 100 h kms" 1 Mpc" 1 
with h — 0.7. In this cosmology, 1' corresponds to 
238kpc/i~ 1 = 341 kpc at the cluster redshift, z = 0.439. 
We use the standard notation Ma = M 3 d(< ta) to de- 
note the total mass enclosed within a sphere of radius 
ta, within which the mean interior density is A times 
the critical mass density at the cluster redshift. We refer 
all our virial quantities to an overdensity A of A v i r ~ 132 
based on the spherical collapse model (Appendix A of 
Kitayama & Suto 1996). 32 All quoted errors are 68.3% 
confidence limits (CL) unless otherwise stated. The ref- 
erence sky position is the center of the brightest cluster 
galaxy (BCG) of Zitrin et al. (2012), R.A. = 12:06:12.15, 
Deck = -08:48:03.4 (J2000.0). 

2. BASIC THEORY OF GALAXY CLUSTER WEAK 
LENSING 

The central quantity of interest in this work is the 
convergence of gravitational lensing, k(0) = £(0)/£ cr it, 
which is the surface mass density projected on to the lens 
plane, £(0), in units of the critical surface mass density 
for lensing, 



AnGDi 



Dis{z s ) 
D s (z s ) 



(1) 



Here D s , Di, and Di s are the proper angular diameter 
distances from the observer to the source, from the ob- 
server to the lens, and from the lens to the source, respec- 
tively; /3 is the angular-diameter distance ratio associated 
with the population of background sources. 

The lens distortion and magnification of images are 
described by the Jacobian matrix A a p (a,/3 = 1,2) of 
the lens mapping, which can be decomposed as A a p = 
(1 — K)S a p — T a p, where 6 a p is Kronecker's delta, and 
T a p is the trace-free, symmetric shear matrix, 



_ ( +71 72 
72 -7i 



(2) 



with the components of complex gravitational shear 
with spin-2 nature (under coordinate rotations; sec 
Bartelmann & Schneider 2001; Okura et al. 2008), 7 = 



7i + *72 = I7|e 
each other by 



The k and 7 fields are related to 



Ak(0) = d a d fi T a f ) (e). 



(3) 



The Green's function for the two-dimensional (2D) Pois- 
son equation is A _1 (0,0') = In |0 - 0'\/(2n), so that 
Equation (3) can be readily solved (Kaiser & Squires 
1993). 

In the strict weak-lensing limit (k, I7I -C 1), r Q ^ in- 
duces a small quadrupolc distortion of the background 
image, which can be measured from observable clliptic- 
ities of background galaxy images (Kaiser et al. 1995). 
In general, the observable quantity for quadrupole weak 

32 A v ; r ft: 134 using the fitting formula given by 
Bryan & Norman (1998). 
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lcnsing is not 7 but the reduced gravitational shear, 

1 

9 = 9i + i92 = (4) 

1 — K 

in the subcritical regime where dct.4 > (or 1/g* in the 
negative parity region with det.4 < 0). 

Given an arbitrary circular loop of radius 9 on the sky, 
the tangential shear 7+ (0) averaged around the loop sat- 
isfies the following identity (e.g., Kaiser 1995): 

1+ {9) = k{<9)-k(9), (5) 

where n{9) is the azimuthal average of k around the loop, 
and 7t(< 9) is the average convergence within the loop. 
Hence, a constant mass sheet cannot be constrained us- 
ing the shear information alone, known as the mass-sheet 
degeneracy (e.g., Bartelmann & Schneider 2001). 

This inherent degeneracy can be unambiguously bro- 
ken by measuring the magnification effects, which pro- 
vide complementary and independent constraints on the 
cluster mass distribution (Umetsu et al. 2011b). The 
magnification is given by the inverse Jacobian determi- 
nant, 



PV ; |deU(0)| |(1- K ) 2 -| 7 | 2 |- 

The magnification ji{0) can influence the observed sur- 
face density n^{6) of background sources, expand- 
ing the area of sky, and enhancing the observed 
flux of background sources (e.g., Broadhurst et al. 
1995; Umetsu & Broadhurst 2008; Van Waerbeke et al. 
2010; Rozo k. Schmidt 2010; Umetsu et al. 2011b; 
Hildebrandt et al. 2011; Ford et al. 2011). The former 
effect reduces the effective observing area in the source 
plane, decreasing the number of sources per solid angle; 
on the other hand, the latter effect amplifies the flux of 
background sources, thereby increasing the number of 
sources above the limiting flux. The net effect is known 
as magnification bias, and depends on the intrinsic slope 
of the luminosity function of background sources as: 

n^B) ^n Qf i(e) 2 - 5s -\ (7) 

where uq = dNo(< m cut )/dfl is the unlensed mean num- 
ber density of background sources for a given magni- 
tude cutoff TOcut, approximated locally as a power- law 
cut with slope, s = d\og 10 N (< m)/dm > 0. In the 
strict weak-lensing limit, n^/n^ — 1 ps (5s — 2)k. For a 
maximally-depleted population of galaxies with s = 0, 
n^/riQ = /.t -1 pa 1 — 2k in this limit. 

Alternatively, the mass-sheet degree of freedom can be 
determined such that the mean S averaged over the out- 
ermost cluster region vanishes, if a sufficiently wide sky 
coverage is available. 33 

3. SUBARU DATA AND ANALYSIS 

In this section we present a technical description of our 
weak-lensing analysis of MACS1206 based on deep Sub- 
aru multi-color images. The data reduction and the pho- 
tometry procedure are summarized in Section 3.1. The 

33 Or, one may constrain the constant such that the enclosed 
mass within a certain aperture is consistent with cluster mass es- 
timates from some other observations (e.g., Umetsu & Futamasc 
2000). 



details of our weak-lensing shape analysis are given in 
Section 3.2. Our shear calibration strategy is described 
in Section 3.3. Details of the sample selection and lens- 
ing depth estimation are given in Sections 3.4 and 3.5, 
respectively. 

3.1. Subaru Data and Photometry 

We analyze deep BVRJ c z' images of MACS1206 ob- 
served with the wide-field camera Suprime-Cam (34' x 
27'; Miyazaki et al. 2002) at the prime focus of the 8.3-m 
Subaru telescope. The observations are available in the 
Subaru archive, SMOKA. 34 The seeing FWHM in the 
co-added mosaic image is 1.01" in B (2.4ks), 0.95" in 
V (2.2ks), 0.78" in R c (2.9 ks), 0.71" in I c (3.6ks), and 
0.58" in z' (1.6ks) with 0.20" pixel" 1 , covering a field of 
approximately 36' x 34'. The limiting magnitudes are ob- 
tained asB = 26.5, V = 26.5, R c = 26.2, I c = 26.0, and 
z' = 25.0 mag for a 3er limiting detection within 2" di- 
ameter aperture. The observation details of MACS1206 
are listed in Table 2. Figure 1 shows Subaru BVR c I c z' 
composite color images of the cluster field, produced au- 
tomatically using the publicly available Trilogy software 
(Coe et al. 2012). 35 

Standard reduction steps were performed using the 
mscred task in IRAF. 36 We closely follow the data re- 
duction procedure outlined in Nonino et al. (2009) to cre- 
ate a co-added mosaic of Subaru Suprime-Cam images, 
incorporating additional reduction steps, such as auto- 
mated masking of bleeding of bright saturated stars. 

To obtain an accurate astrometric solution for Sub- 
aru observations, we retrieved processed MegaCam griz 
images from the CFHT archive, 37 and used MegaCam 
r data (Filter Number: 9601) as a wide-field reference 
image. A source catalog was created from the co-added 
MegaCam r image, using the 2MASS catalog 38 as an 
external reference catalog. The extracted r catalog has 
been used as a reference for the SCAMP software (Bertin 
2006) to derive an astrometric solution for the Suprime- 
Cam images. 

The photometric zero points for the co-added Suprime- 
Cam images were bootstrapped from a suitable set of 
reference stars identified in common with the calibrated 
MegaCam data. These zero-points were refined in two 
independent ways: firstly by comparing cluster elliptical- 
type galaxies with the HST/ ACS images; subsequently 
by fitting SED (spectral energy distribution) templates 
with the BPZ code (Bayesian photometric redshift es- 
timation, Bcm'tez 2000; Benftez et al. 2004) to Subaru 
photometry of 1163 galaxies having measured spectro- 
scopic redshifts from VLT/VIMOS (P. Rosati et al., in 

34 http://smoka.nao.ac.jp 

35 http:/ /www. stsci.edu/~dcoc/trilogy/ 

36 IRAF is distributed by the National Optical Astronomy Ob- 
servatories, which is operated by the Association of Universities for 
Research in Astronomy, Inc., under cooperative agreement with the 
National Science Foundation 

37 This research used facilities of the Canadian Astronomy Data 
Centre operated by the National Research Council of Canada with 
the support of the Canadian Space Agency 

38 This publication makes use of data products from the Two 
Micron All Sky Survey, which is a joint project of the Univer- 
sity of Massachusetts and the Infrared Processing and Analysis 
Center / California Institute of Technology, funded by the National 
Aeronautics and Space Administration and the National Science 
Foundation 
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Figure 1. Subaru BVR c I c z' composite color images centered on the galaxy cluster MACS1206 (z = 0.439), overlaid with mass contours 
from our joint strong-and-wcak lensing analysis (SaWLens) of HST and Subaru observations. The image size in the left panel is 24' X 24' 
covering a projected area of 5.7 X 5.7Mpch -2 at the cluster redshift. In the left and right panels, the lowest contour levels are re = 0.12 
and 0.15, with increments of Are = 0.09 and 0.07, respectively. The right panel is a zoom-in- view of the boxed region of the left panel, with 
a side length of 8' (1.9 Mpc h -1 ). North is top and east is left. 

preparation). This leads to a final photometric accu- 
racy of ~ 0.01 mag in all five passbands (see also Section 
3.5). Five-band BVR c I c z' photometry catalog was then 0M 
measured using SExtractor (Bertin & Arnouts 1996) in 
point-spread-function (PSF) matched images created by 
ColorPro (Coe et al. 2006), where a combination of all *<& 

five bands was used as a deep detection image. The stel- 
lar PSFs were measured from a combination of 100 stars 
per band and modeled using IRAF routines. -o.oa 

For the weak- lensing shape analysis (Section 3.2), we 
use the J c -band data taken in 2009 January, which have -0.04 o o.M -0.04 o o.M 

the best image quality in our data-sets (in terms of the *i *i 

stability and coherence of the PSF anisotropy pattern, 
taken in fairly good seeing conditions). Two separate 
co-added 7 c -band images, each with a total exposure 
time of 1.1 ks, were produced based on the imaging ob- 
tained at two different camera orientations separated by 
90 degrees, in order not to degrade the shape measure- 
ment quality. 



3.2. Subaru Weak Lensing Shape Analysis 

For shape measurements of background galaxies, we 
use our weak-lensing analysis pipeline based on the IM- 
CAT package (Kaiser et al. 1995, KSB hereafter), in- 
corporating modifications and improvements outlined in 
Umetsu et al. (2010). Our KSB+ implementation has 
been applied extensively to Subaru cluster observations 
(e.g., Broadhurst et al. 2005a, 2008; Umetsu et al. 2007; 
Umetsu & Broadhurst 2008; Okabe & Umetsu 2008b; 
Umetsu et al. 2009, 2010, 2011b,a; Medezinski et al. 
2010, 2011; Zitrin et al. 2011c; Coe et al. 2012). 

We measure components of the complex image ellip- 
ticity, e a = {Qu - Q22, Q12} /(Qu + Q22), from the 
weighted quadrupole moments of the surface brightness 
1(6) of individual objects, 



Q a0 = d 2 ew(0)e a p i(e) (0,0 = 1,2) (8) 



Figure 2. Stellar ellipticity distributions before and after the 
PSF anisotropy correction for Subaru/Suprimc-Cam J c -band data 
taken with camera orientations of PA = 0° (Orientation 1; red) 
and PA = 90° (Orientation 2; black). The left panel shows the 
raw ellipticity components (ej, e^) of stellar objects, and the right 
panel shows the residual ellipticity components (Se*, Se^) after the 
PSF anisotropy correction. 

where W(9) is a Gaussian window function matched to 
the size (r g ) of the object, and the weighted object cen- 
troid is chosen as the coordinate origin, which is itera- 
tively refined to accurately measure the object shapes. 

Next, we correct observed ellipticities e a for the PSF 
anisotropy using a sample of stars in the field as ref- 
erences. We select bright (18 I c ~ 22), unsaturated 
stellar objects identified in a branch of the object half- 
light radius (r/j) versus I c diagram, and measure the 
PSF anisotropy kernel of the KSB algorithm as a func- 
tion of the object size r g . Figure 2 shows the distri- 
butions of stellar ellipticity components (e* ) before and 
after the PSF anisotropy correction. From the rest of 
the object catalog, we select as a weak-lensing galaxy 
sample those objects with v > 10, Th > tT + 1.5cr(rf ), 
and r g > modc(r*), where v is the KSB detection sig- 
nificance, r£ and cr(r^) are median and rms dispersion 
values of stellar sizes rj*. The anisotropy corrected ellip- 
ticities e' a arc then corrected for the isotropic smearing 
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effect as g a = e' a /P g . 
For each galaxy we assign the statistical weight, 

W (k) = 2 -T 2 ' ( 9 ) 

fl(fc) g 

where c 2 ( fc ) is the variance for the reduced shear estimate 
of the fcth galaxy computed from 50 neighbors identified 
in the r g -I c plane, and a 2 g is the softening constant vari- 
ance (e.g., Hamana ct al. 2003; Umetsu & Broadhurst 
2008; Oguri et al. 2009; Okabe et al. 2010). This weight- 
ing scheme is essential to down-weight faint and small 
objects which have noisy shape measurements (see Fig- 
ure 4 of Umetsu et al. 2010). We choose a g = 0.4, which 
is a typical value of the mean rms a g over the back- 
ground sample (see Table 3; Umetsu & Broadhurst 2008; 
Umetsu et al. 2009; Okabe et al. 2010). 
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3.3. Shear Calibration 

We follow the shear calibration strategy of 
Umetsu et al. (2010) to improve the precision in 
shear recovery. This is motivated by the general ten- 
dency of KSB+ to systematically underestimate the 
shear signal in the presence of measurement noise (see 
Umetsu et al. 2010; Okura & Futamase 2011). 

First, we select as a sample of shear calibrators those 
galaxies with v > v c and P g > 0. Here we take v c = 20. 
Note that the shear calibrator sample is a subset of the 
target galaxy sample. Second, we divide the calibrator 
r g -I c plane into a grid of 2 x 10 cells each containing 
approximately equal numbers of calibrators, and com- 
pute a median value of P g at each cell. Then, each ob- 
ject in the target sample is matched to the nearest point 
on the (r g ,I c ) calibration grid to obtain a filtered mea- 
surement, (P g ). Finally, we use the calibrated estimator 
9a = e' a /(Pg) f° r the reduced shear. 

We have analyzed the two I c mosaic images separately 
to construct a composite galaxy shape catalog, by prop- 
erly weighting and combining the calibrated distortion 
measurements (g a ) for galaxies in the overlapping region. 

We have tested our analysis pipeline using simu- 
lated Subaru Suprime-Cam images (see Section 3.2 of 
Oguri et al. 2012; Massey et al. 2007). We find that we 
can recover the weak-lensing signal with good precision, 
typically, \m\ < 5% of the shear calibration bias, where 
the range of m- values shows a modest dependence of cal- 
ibration accuracy on seeing conditions and PSF proper- 
ties, and c ~ 10~ 3 of the residual shear offset, which 
is about one order of magnitude smaller than the typ- 
ical distortion signal in cluster outskirts (\g\ ~ 10 -2 ). 
This level of performance is comparable to other similarly 
well-tested methods (Heymans ct al. 2006; Massey et al. 
2007). 

3.4. Sample Selection 

A careful background selection is critical for a 
wcak-lcnsing analysis so that unlensed cluster mem- 
bers and foreground galaxies do not dilute the true 
lensing signal of the background (Broadhurst et al. 
2005a; Medezinski et al. 2007; Umetsu & Broadhurst 
2008; Medezinski et al. 2010). This dilution effect is sim- 
ply to reduce the strength of the lcnsing signal when av- 
eraged over a local ensemble of galaxies (by a factor of 



Figure 3. Blue and red background galaxies are selected for 
weak-lensing analysis (lower left blue dashed and right red dot- 
dashed regions, respectively) on the basis of Subaru BR c z' color- 
color-magnitude selection. All galaxies with z 1 < 24.6 mag (cyan) 
are shown in the diagram. At small radius, the cluster overdensity 
is identified as the green outlined region, defining our green sample 
comprising mostly the red sequence of the cluster and a blue trail 
of later type cluster members. The background samples are well 
isolated from the green region and satisfy other criteria as discussed 
in Section 3.4. Our background selection successfully excludes all 
spectroscopically-confirmed cluster members (black) found within 
the projected cluster virial radius (r v i r ks 1.6 Mpc h~ 1 ). The 
cluster members are determined from the ongoing survey with 
VLT/VIMOS (Rosati ct al., in prep.), using the algorithm of 
Mamon et al. (2010) in the dynamical analysis which will be pre- 
sented in a forthcoming paper (Biviano et al. in prep.). 
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Figure 4. Surface number density profiles n(9) of Subaru BR c z'- 
selected galaxies used for the weak-lensing shape analysis. The 
results arc shown for our red (triangles), blue (circles), and green 
(crosses) samples. See also Figure 9. 

2-5 at R £ 400 kpc hr 1 ; see Figure 1 of Broadhurst et al. 
2005a), particularly at small cluster radius where the 
cluster is relatively dense, in proportion to the fraction 
of unlensed galaxies whose orientations are randomly dis- 
tributed. 

We use the background selection method of 
Medezinski ct al. (2010) to define undiluted sam- 
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Figure 5. Azimuthally-avcraged radial profiles of the tangential 
reduced shear <;+ (upper panel) and the 45° rotated (x) component 
g x (lower panel) for our Subaru red (triangles), blue (circles), green 
(crosses), and blue+red (squares) galaxy samples shown in Figure 
4. The error bars represent 68.3% confidence intervals estimated by 
bootstrap resampling techniques. The symbols for the red and blue 
samples arc horizontally shifted for visual clarity. For a consistency 
check, we compare our Subaru results with CFHT/Megacam data 
based on our grz-selected background sample (gray area). The 
profile for the green sample is consistent with a null signal at all 
radii, while this population is strongly clustered at small radius 
(Figures 4), indicating that the green galaxies mostly consist of 
cluster member galaxies. For all of the samples, the X -component 
is consistent with a null signal detection well within 2a at all radii, 
indicating the reliability of our distortion analysis. 



pies of background galaxies, which relies on empirical 
correlations for galaxies in color-color-magnitude space 
derived from the deep Subaru photometry, by reference 
to evolutionary tracks of galaxies (for details, see 
Medezinski et al. 2010; Umetsu et al. 2010) as well as to 
the deep photometric-redshift survey in the COSMOS 
field (Ilbert et al. 2009). 

For MACS1206, we have a wide wavelength coverage 
(BV R c I c z') of Subaru Suprime-Cam. We therefore make 
use of the (B—R c ) vs. (R c — z') color-color (CC) diagram 
to carefully select two distinct background populations 
which encompass the red and blue branches of galaxies. 
We limit the data to z' = 24.6 mag in the reddest band, 
corresponding approximately to a 5c limiting magnitude 
within 2" diameter aperture. Beyond this limit incom- 
pleteness creeps into the bluer bands, complicating color 
measurements, in particular of red galaxies. 

To do this, we first identify in CC space an overdensity 
of galaxies with small projected distance < 3' ( ^ 1 Mpc 
at zi = 0.439) from the cluster center. Then, all galax- 
ies within this distinctive region define the green sample 
(see the green outlined region in Figure 3), comprising 
mostly the red sequence of the cluster and a blue trail 
of later type cluster members (Medezinski et al. 2010; 
Umetsu et al. 2010), showing a number density profile 
that is steeply rising toward the center (Figure 4, green 
crosses). The weak-lensing signal for this population is 
found to be consistent with zero at all radii (Figure 5, 
green crosses), indicating the reliability of our procedure. 
For this population of galaxies, we find a mean photomet- 
ric rcdshift of (z p hot) ~ 0.44 (see Section 3.5), consistent 
with the cluster rcdshift. Importantly, the green sample 
marks the region that contains a majority of unlensed 



galaxies, relative to which we select our background sam- 
ples, as summarized below. 

For the background samples, we define conservative 
color limits, where no evidence of dilution of the weak- 
lensing signal is visible, to safely avoid contamination 
by unlensed cluster members and foreground galaxies. 
The color boundaries for our blue and red background 
samples are shown in Figure 3. For the blue and red 
samples, we find a consistent, clearly rising weak-lensing 
signal all the way to the center of the cluster, as shown 
in Figure 5. 

For validation purposes, we compare in CC space 
our color samples with a spectroscopic sample of 
cluster galaxies in MACS1206. Figure 3 shows 
that the background selection procedure established in 
our earlier work (Medezinski et al. 2010; Umetsu et al. 
2010; Medezinski et al. 2011) successfully excludes 
all spectroscopically-confirmcd cluster members found 
within the projected cluster virial radius (r V j r rj 
1.6 Mpc ft -1 ; see Section 6). The cluster members are 
determined from the ongoing survey with VLT/VIMOS, 
part of the VLT-CLASH Large Programme 186.A-0798 
(P. Rosati et al., in preparation), using the algorithm 
of Mamon et al. (2010) in the dynamical analysis which 
will be presented in a forthcoming paper (A. Biviano et 
al. in preparation). We find about 70% of the clus- 
ter members overlap with our CC-selected green galax- 
ies; the rest are cluster members with bluer colors. We 
note there is a statistically inevitable fraction of interlop- 
ers even in the dynamically-selected cluster membership 
as discussed in Wojtak et al. (2007, their Table 1) and 
Mamon et al. (2010, their Figure 13). 

As a further consistency check, we also plot in Fig- 
ure 4 the galaxy surface number density BjS ci function 
of radius, n(8), for the blue and red samples. As can 
be seen, no clustering is observed toward the center for 
the background samples, which demonstrate that there 
is no significant contamination by cluster members in the 
samples. The red sample reveals a systematic decrease in 
their projected number density toward the cluster cen- 
ter, caused by the lensing magnification effect (Section 
2). A more quantitative magnification analysis is given 
in Section 5.2.2. 

To summarize, our CC-selection criteria yielded a total 
of N = 13252, 1638, and 4570 galaxies, for the red, green, 
and blue photometry samples, respectively (Table 3). For 
our weak-lensing distortion analysis, we have a subset of 
8969 and 4154 galaxies in the red and blue samples (with 
usable I c shape measurements), respectively (Table 4). 

3.5. Depth Estimation 

The lensing signal depends on the source redshift z s 
through the distance ratio /3(z s ) = Di s /D s . We thus 
need to estimate and correct for the respective depths 
of the different galaxy samples, when converting the 
observed lensing signal into physical mass units. 

For this we used BPZ (Section 3.1) to measure pho- 
tometric rcdshifts (photo-zs) z p i lot for our deep Sub- 
aru BVR c I c z' photometry (Section 3.1). BPZ employs 
a Baycsian inference where the redshift likelihood is 
weighted by a prior probability, which yields the proba- 
bility density P(z,T\m) of a galaxy with apparent mag- 
nitude m of having certain redshift z and spectral type 
T. In this work we used a new library (N. Benitez 2012, 
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in preparation) composed of 10 SED templates originally 
from PEGASE (Fioc & Rocca-Volmerange 1997) but re- 
calibrated using the FIREWORKS photometry and spec- 
troscopic redshifts from Wuyts et al. (2008) to optimize 
its performance. This library includes five templates for 
elliptical galaxies, two for spiral galaxies, and three for 
starburst galaxies. In our depth estimation we utilize 
BPZ's ODDS parameter, which measures the amount of 
probability enclosed within a certain interval Az centered 
on the primary peak of the redshift probability density 
function (PDF), serving as a useful measure to quan- 
tify the reliability of photo- z estimates (Bemtez 2000). 39 
We used our VLT/VIMOS sample of 1163 galaxies with 
spectroscopic redshifts z spcc ( £ 1.5) to assess the perfor- 
mance of our photo- z estimation. From the whole sam- 
ple, we find an rms scatter of a(8 z ) « 0.027 in the frac- 
tional error S z = (z phot - z spoc )/(l + z spoc ), with a small 
mean offset n(b~ z ) = —0.0021 and a 5<j outlier fraction 
of « 5.5%. Using a subsamplc of ~ 510 galaxies with 
0.3 < z spcc < 0.5, we find a{8 z ) 0.031 with w 1.5% of 
outliers. 

For a consistency check, we also make use of the COS- 
MOS catalog (Ilbert et al. 2009) with robust photometry 
and photo- z measurements for the majority of galaxies 
with i' < 25 mag. For each sample, we apply the same 
CC selection to the COSMOS photometry, and obtain 
the redshift distribution N(z) of field galaxies. 

For each background population, we calculate weighted 
moments of the distance ratio j3 as 

_ fdz W (z)N(z)(3 n (z) 
{l ' Jdzw{z)N(z) ' 1 j 

where w(z) is a weight factor; w is taken to be the 
Bayesian ODDS parameter for the BPZ method, and 
w = 1 otherwise. The sample mean redshift (z s ) is de- 
fined similarly to Equation (10). The first moment (/3) 
represents the mean lensing depth. 40 It is useful to define 
the effective single-plane source redshift, z Sje ff, such that 
(Umetsu & Broadhurst 2008; Umetsu et al! 2009, 2010) 

/?(z s>eff ) = (/?>. (11) 

In Table 4 we summarize the mean depths (j3) and 
the effective source redshifts z Si0 ff for our background 
samples. For each background sample, we obtained con- 
sistent mean-depth estimates ((3) (within 2%) using the 
BPZ- and COSMOS-based methods. In the present 
work, we adopt a conservative uncertainty of 5% in the 
mean depth for the combined blue and red sample of 
background galaxies, (/3(back)) = 0.54±0.03, which cor- 
responds to z SjC ff = 1.15 ± 0.1. We marginalize over this 
uncertainty when fitting parameterized mass models to 
our weak-lcnsing data. 

4. CLUSTER STRONG LENSING ANALYSIS 

39 In the present work, we set Az = 2 X 0.03(1 + z p hot)> which 
is approximately twice the width (a) of the error distribution. 

40 In general, a wide spread of the redshift distribution of back- 
ground galaxies, in conjunction with the single-plane approxima- 
tion, may lead to an overestimate of the gravitational shear in the 
nonlinear regime (Hoekstra et al. 2000). Thanks to the deep Sub- 
aru photometry, we found that this bias in the observed reduced 
shear is approximately Ag/g Kb ((/3 2 )/(/3) 2 — l)re ~ 0.06k to the 
first order of re. See § 3.4 of Umetsu et al. (2010) for details. 



For a massive cluster, the strong- and weak-lcnsing 
regimes contribute quite similar logarithmic coverage of 
the radial mass profile. It is therefore crucial to include 
the central strong-lcnsing information in a cluster lensing 
analysis (e.g., Umetsu et al. 201 lb, a). 

Here we perform several complementary strong-lcnsing 
analyses using a wide variety of modeling methods, 
namely the Zitrin et al. (2009) method, Lenstool 
(Kneib et al. 1996; Julio et al. 2007), LensPerfect 
(Coe et al. 2010, 2012), Pixelens (Saha & Williams 
2004; Grillo et al. 2010), and a joint strong-and-weak 
lensing reconstruction method of Merten et al. (2009, 
2011b) (hereafter, SaWLens). All analyses here use the 
positions and redshifts of multiply-lensed images identi- 
fied by Zitrin et al. (2012). 

Lens reconstruction methods are broadly classified into 
parametric and non-parametric: In the former approach, 
the total mass distribution of the deflector is described 
in terms of a set of theoretically (and/or observationally) 
motivated models, each specified by a particular func- 
tional form characterized by a small number of free pa- 
rameters. This involves, to some extent, the assignment 
of halos to visible galaxies assuming light approximately 
traces mass, while the latter does not except for certain 
priors on the mass distribution. 41 Among the methods 
used in the present work, the Zitrin et al. (2009) method 
and Lenstool arc parametric; LensPerfect, Pixe- 
lens, and SaWLens are non-parametric. 

For this work, we primarily use the detailed strong lens 
modeling of Zitrin et al. (2012) based on deep CLASH 
imaging and VLT/VIMOS spectroscopy, as summarized 
in Section 4.1. The cluster miscentering effects are dis- 
cussed in Section 4.2. In Section 4.3 we introduce and 
apply a technique to self-calibrate the bin-bin covariance 
matrix of the central radial mass profile derived from the 
reanalysis of Zitrin et al. (2012). In Section 4.4 we per- 
form several semi-independent strong-lensing analyses on 
the MACS1206 HST images, utilizing various modeling 
methods, in order to verify the identifications of the mul- 
tiple images and to independently assess the level of in- 
herent systematic uncertainties in our analyses. 

4.1. Primary Strong Lensing Model 

Here we briefly summarize our well-tested approach 
to strong-lens modeling, developed by Broadhurst et al. 
(2005b) and optimized further by Zitrin et al. (2009), 
which has previously uncovered large numbers of 
multiply-lensed galaxies in HST images of many clusters 
(e.g., Broadhurst et al. 2005b; Zitrin et al. 2009, 2010, 
2011a,b,c). In the present work, we use a new Markov 
Chain Monte Carlo (MCMC) implementation of the 
Zitrin et al. (2009) method, where also the BCG mass 
is allowed to vary. 42 

Our flexible mass model consists of four components, 
namely the BCG, cluster galaxies, a smooth DM halo, 
and the overall matter ellipticity (corresponding to a 
coherent external shear; for details, see Zitrin et al. 
2009), described by 7 free parameters in total. 43 The 

41 The latter is often based on the assumption that the lens pro- 
files and/or distributions can be well approximated by a pixelatcd 
mass distribution (e.g., Pixelens and SaWLens). 

42 Our very preliminary MCMC results were presented in Figure 
4 of Zitrin et al. (2012). 

43 The Zitrin et al. (2009) method employs grid-based 
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basic assumption adopted is that cluster galaxy light 
approximately traces the DM; the latter is modeled 
as a smoothed version of the former (see, for details, 
Zitrin et al. 2009). This approach to strong lensing is 
sufficient to accurately predict the locations and internal 
structure of multiple images, since in practice the number 
of multiple images uncovered readily exceeds the number 
of free parameters, so that the fit is fully constrained. 

Zitrin et al. (2012) identified 47 new multiple images of 
12 distant sources (including 3 candidate systems; Sys- 
tems 9-11 therein), in addition to the known giant arc 
system at z s — 1.03 (Ebeling et al. 2009), bringing the 
total known for this cluster to 50 multiply-lensed im- 
ages of 13 sources, spanning a wide redshift range of 
1 ^ z s ^ 5.5, spread fairly evenly over the central region, 
3" < < 1'. Zitrin et al. (2012) used the position and 
redshift of 32 secure multiple images of 9 systems to con- 
strain the mass model. Following Zitrin et al. (2012), we 
adopt an image positional error of 2" (w 1.4" in each 
dimension), which is a typical value in the presence of 
uncorrelated large scale structure (LSS) along the line 
of sight (for details, see Zitrin et al. 2012; Host 2012; 
Julio et al. 2010). Including the BCG mass as an addi- 
tional free parameter, we find here an acceptable fit with 
the minimized % 2 value (Xmin) OI 22.8 for 39 degrees of 
freedom (dof), with an image-plane reproduction error 
of 1.76". The new MCMC results are in good agreement 
with the results of Zitrin et al. (2012), as shown here in 
Figure 6, with only some minor differences at the inner- 
most radii £ 2" (~ 8kpcft. _1 ) dominated by the BCG 
(see Newman et al. 2009). The detailed central mass 
map reveals a fairly elliptical outer critical curve (see 
Figure 1 of Zitrin et al. 2012). For a source at z s = 2.54, 
the outer critical curve encloses an area with an effective 
Einstein radius of 6*Ein = 28" ± 3"; for the lower- redshift 
system with z s = 1.03, the effective Einstein radius of 
the critical area is 0Ein = 17" ± 3" (Table 1). 

4.2. Cluster Miscentering Effects 

To obtain meaningful radial profiles, one must care- 
fully define the cluster center. It is often assumed that 
the cluster mass centroid coincides with the BCG po- 
sition, whereas BCGs can be offset from the mass cen- 
troids of the corresponding DM halos (Johnston et al. 
2007; Oguri & Takada 2011; Umetsu et al. 2011b,a). 

Here we utilize our detailed mass model of Zitrin et al. 
(2012), which allows us to locate the peak position of the 
smooth DM component, providing an independent mass 
centroid determination (e.g., Umetsu et al. 2010, 2011b). 
In this method, we approximate the large-scale distribu- 
tion of cluster mass by assigning a power-law mass profile 
to each cluster galaxy, the sum of which is then smoothed 
to represent the DM distribution. The success of this 
simple model in describing the projected mass distribu- 
tions of lensing clusters, as well as identifying many sets 
of multiply-lensed images, assures us that the effective 
DM center can be determined using multiple images as 
well as the distribution of cluster member galaxies. In 
this context, the DM peak location is primarily sensi- 
tive to the degree of smoothing (S) and the index of the 

maximum-likelihood parameter estimation in the 6-dimensional pa- 
rameter space, where the seventh parameter included in the present 
work is the BCG mass. 
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Figure 6. Surface mass density profile n derived from our Subaru 
weak-lcnsing and Hubble strong-lensing measurements. The red 
circles represent our full weak-lcnsing constraints from joint shear 
and magnification measurements (Figure 9), consistent with the 
purely shear-based results (squares) and the SaWLens results (or- 
ange line with error bars), all showing a shallow radial trend with 
a nearly isothermal logarithmic density slope, dlnS/dlni?, ~ — 1. 
For weak lensing, the innermost bin represents the average con- 
vergence 7t(< m i n ) interior to the inner radial boundary of the 
weak-lensing data (0.8' < 9 < 16'), m i n = 0.8', which is about 
twice the Einstein radius for a distant background source at z s ~ 2 
(see Tables 3 and 4), and hence sufficiently large for our back- 
ground galaxies at an effective source redshift of z s e jf = 1.15 ±0.1. 
The triangles show the NE-SW mass profile excluding the large 
scale structure extending along the NW-SE direction (see Figure 
8), derived from a two-dimensional mass reconstruction using both 
shear and magnification data, in good agreement with the stan- 
dard NFW form (gray area). The black solid line is the best-fit 
model of Zitrin et al. (2012) based on the grid-based maximum 
likelihood parameter estimation. The small blue circles with error 
bars represent our primary strong-lens constraints on the binned 
mass profile derived from an MCMC implementation of Zitrin et al. 
(2012). The errors are based on the self-calibrated covariance ma- 
trix (only every other point is shown here; Section 4.3). Our mass 
profile results from several weak and strong lensing methods all 
agree in the regions of overlap within their corresponding uncer- 
tainties. For the sake of clarity, the Pixelens and Lenstool results 
are shown without error bars. The bottom panel shows the respec- 
tive deviations Ak (in units of the error a) from the best-fit NFW 
model. The projected mass profile averaged over all azimuthal an- 
gles (squares, circles) shows a systematic excess at large radii with 
R £ IMpcft- 1 (6» > 4'). 

power law (q) of Zitrin et al. (2009). 

We find only a small offset of ~ 1", or a projected off- 
set distance of d Q s = 4kpc/i _1 at zi = 0.439, between 
the BCG and the DM peak of mass, well within the un- 
certainties. The BCG position also coincides well with 
the peak of X-ray emission within 2" in projection (Table 
1). This level of cluster centering offset is fairly small as 
compared to those found in other high-mass clusters, say 
doff ps 20kpc/i~ 1 in RXJ1347-11 (Umetsu et al. 2011b), 
often implied by other massive bright galaxies in the 
vicinity of the BCG. In the present work, we thus adopt 
the BCG position as the cluster center, and limit our 
analysis to radii greater than 4" (« 16kpc/i _1 ), which 
is approximately the location of the innermost strong- 
lensing constraint (see Section 4.1) and sufficiently large 
to avoid the BCG contribution. This inner radial limit 
corresponds roughly to 4d ff(> 2<i ff), beyond which 
smoothing from the cluster miscentering effects on the 
£ profile is sufficiently negligible (Johnston et al. 2007; 
Umetsu et al. 2011a; Sereno & Zitrin 2012). 
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4.3. Self- Calibration of the Covariance Matrix 

The MCMC approach allows for a full parameter-space 
extraction of the underlying lensing signal. We construct 
from MCMC samples a central mass profile Ki and its co- 
variance matrix Cij in linearly-spaced radial bins, span- 
ning from 6 = 1" to the limit of our ACS data, 6 ~ 100". 
Note, multiple image constraints are available out to a 
radius of w V (Section 4.1), so that the mass model be- 
yond this radius is constrained essentially by the light 
distribution of cluster member galaxies, and hence the 
constraints there are driven by the prior. We find that 
the mass profile is positively correlated from bin to bin, 
especially at radii beyond #Ein ~ 28" (z s = 2.5). Accord- 
ingly, the C matrix is nearly singular, with very small 
eigenvalues associated with large-scale modes where the 
constraints are weaker, leading to underestimated diag- 
onal errors at £ 6> E i n ~ 28" (z s = 2.5). 

Here, we use a regularization technique with a single 
degree of freedom to calibrate the C matrix and obtain 
conservative errors for strong lcnsing, accounting for pos- 
sible systematic errors introduced by the prior assump- 
tions in the modeling. We first perform an eigenvalue de- 
composition as C = [/AC/*, where A is a diagonal matrix 
of eigenvalues and U is a unitary matrix of eigenvectors. 
Then, we determine our regularization constant, the min- 
imum eigenvalue A m i n , by conservatively requiring that 
the outermost k value, K m i n = «(100") ~ 0.22, is consis- 
tent with a null detection: i.e., A m ; n = K^ in = (0.22) 2 . 
Replacing those less than A m j n by A m j n and restoring the 
C matrix with the regularized A yields the desired, self- 
calibrated C matrix. All points at ^ 1' are then excluded 
from our analysis. We find that a weaker regulariza- 
tion with A m j n = (0.1) 2 only affects the halo parameters 
(M v ; r , c v ir) by less than 4%. 

In Figure 6 we show our strong-lensing constraints on 
the central K profile using the self-calibrated C matrix, 
where the outer radial boundary is conservatively set to 
9 = 53" (w 26> E in at z s = 2; see Zitrin et al. 2012). This 
calibration scheme produces conservative error estimates. 
Overall, the level of correction applied to the C matrix 
increases with increasing radius. We introduce here an 
estimator for the total signal-to-noise ratio (S/N) for de- 
tection, integrated over the radial range considered, and 
quantify the significance of the reconstruction, by the 
following equation (Umetsu & Broadhurst 2008): 

(S/N) 2 = /<C, = k'C- 1 *. (12) 

With the calibrated C matrix, we find a total S/N of 
~ 18 for our strong-lensing k profile in the radial range 
6 < 53". We check that our results are insensitive to the 
choice of radial binning scheme when the self-calibration 
technique is applied. 

4.4. Complementary Strong Lensing Analyses 

We have performed complementary semi-independent 
strong-lensing analyses (Lenstool, LensPerfect, 
Pixelens, SaWLens), using as input the sets (or sub- 
sets) of multiple images identified by Zitrin et al. (2012) 
as well as the same spectroscopic and photometric red- 
shift information. 

In our Lenstool analysis, we parameterize the lens 
mass distribution £(<?) as a multi-component model con- 



sisting of an elliptical NFW potential and truncated ellip- 
tical halos (Kassiola & Kovner 1993) for the 86 brightest 
cluster members. All nine of the secure image systems 
are included as observational constraints. Our best solu- 
tion reproduces all arc systems included and the critical 
lines at z s = 2.54 and 1.03 derived in Zitrin et al. (2012), 
with an image-plane rms of 1.9", very similar to the value 
of ~ 1.8" obtained by Zitrin et al. (2012) and typical to 
parametric mass models for clusters with many multi- 
ple images (Broadhurst et al. 2005b; Halkola et al. 2006; 
Limousin et al. 2007; Zitrin et al. 2009). 

In the Pixelens analysis we model the lens mass dis- 
tribution on a circular grid of 52" radius divided into 18 
pixels. We consider 200 models with decreasing projected 
mass profiles (i.e., oc R~ a with a > 0). We use as 

constraints the spectroscopically-confirmcd Systems 1 to 
4 (Zitrin et al. 2012) of 14 multiple images, spanning the 
range 3.5"to 46"in radius. We check that adding other 
multiple-image systems identified in Zitrin et al. (2012) 
does not significantly affect the Pixelens mass recon- 
struction. 

In the LensPerfect analysis, we assume a prior that 
the projected mass is densest near the center of the 
BCG and decreases outward. Other priors include overall 
smoothness and approximate azimuthal symmetry (for 
details, see Coe et al. 2010). All secure image systems 
are used in this modeling, where including the three can- 
didate systems (#9-11) does not change the results sig- 
nificantly. 

The SaWLens method combines central strong- 
lensing constraints from multiple-image systems with 
weak-lensing distortion constraints in a non-parametric 
manner to reconstruct the underlying lensing potential 
on an adaptively refined mesh. For this cluster we use 
two levels of refinement, providing a 6" pixel resolution 
in the strong-lensing regime covered by CLASH imag- 
ing and a w 22" resolution in the Subaru weak-lensing 
field where the background source galaxies are sparsely 
sampled. The field size for the reconstruction is 25' on a 
side. All image systems except #10 and #11 are in- 
cluded as strong-lensing constraints. The lens distor- 
tion measurements for the blue+red sample are used as 
weak-lensing constraints. The reconstruction errors are 
derived from 1000 bootstrap realizations of the weak- 
lensing background catalog and 1000 samples of the red- 
shift uncertainties in the catalog of strong lensing fea- 
tures. The number of realizations is limited by runtime 
constraints. 

Figure 7 shows and compares the resulting projected 
integrated mass profiles M2d(< 6) derived from our 
comprehensive strong-lensing analyses, along with our 
primary strong-lensing results and model-independent 
Einstein-radius constraints based on Zitrin et al. (2012). 
All these models are broadly consistent with the Ein- 
stein radius constraints. The calibrated error bars 
of Zitrin et al. (2012) are roughly consistent with the 
spread of the semi-independent mass profiles derived 
here. This comparison shows clear consistency among a 
wide variety of analysis methods with different assump- 
tions and systematics, which firmly supports the reliabil- 
ity of our strong-lensing analyses and calibration. 

5. CLUSTER WEAK LENSING ANALYSIS 
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Figure 7. Comparison of projected cumulative mass profiles M2D of MACS1206 obtained with different lensing methods. The red 
shaded area shows our full wcak-lensing constraints (68% CL) derived from a joint Bayesian analysis of Subaru shear and magnification 
measurements (Figure 9) , in good agreement with the shear aperture mass measurements (squares) obtained with a zero density boundary 
condition of E(16' < < 18') = 0. The triangles denote the mass profile using the NE-SW S profile of Figure 6 excluding the NW-SE 
excess regions. The two open rectangles represent model-independent Einstein-radius constraints of 6>Ein = 17" ± 2" (z a = 1.03) and 
0Ein = 28" ± 3" (z a = 2.54). The blue shaded area represents our primary strong-lens model with lcr uncertainty from an MCMC 
implementation of Zitrin et al. (2012), which is broadly consistent with our semi-independent results from a wide variety of four strong-lens 
modeling analyses (Pixelens, Lenstool, LensPerfect, and SaWLens). providing a valuable consistency check. Our independent strong 
and weak lensing profiles are in good agreement in the region of overlap, and together are well described by the standard NFW form (gray 
area), but increasingly exceed it at R IMpc/i -1 out to the limit of our data. The bottom panel shows fractional deviations (AM/M)2r> 
of projected mass profiles with respect to the best-fit NFW model (top, gray), demonstrating the presence of a large scale anisotropy in 
the mass distribution around the cluster. 



This section is devoted to our cluster weak-lensing 
analysis based on the deep multi-color Subaru observa- 
tions. In Section 5.1 we present the projected mass and 
galaxy distributions in and around MACS1206. In Sec- 
tion 5.2 we derive cluster lens distortion and magnifi- 
cation radial profiles from Subaru data. In Section 5.3 
we briefly summarize our Bayesian mass inversion meth- 
ods based on combined lens distortion and magnification 
measurements, and apply to Subaru weak-lensing obser- 
vations of MACS1206. 

5.1. Two- Dimensional Mass Map 

Weak-lensing distortion measurements (g) can be used 
to reconstruct the underlying projected mass density 
field (see Equation 3). Here we use the linear map- 
making method outlined in Section 4.4 of Umetsu et al. 
(2009) to derive the projected mass distribution from the 
Subaru distortion data presented in Section 3. 

In the left panel of Figure 8, we show the £(0) 
field in the central 24' x 24' region, reconstructed from 
the blue+red sample (Section 3.4), where for visualiza- 
tion purposes the mass map is smoothed with a Gaus- 
sian with 1.5' FWHM. A prominent mass peak is vis- 



ible in the cluster center. This first maximum in the 
mass map is detected at a significance level of 9.5<r, 
and coincides well with the optical/X-ray cluster center 
within the statistical uncertainty: AR.A. = 7.0" ± 7.2", 
ADecl. = -1.4" ± 7.6", where AR.A. and ADecl. arc 
right-ascension and declination offsets, respectively, from 
the BCG center. 

Also compared in Figure 8 are member galaxy distribu- 
tions in the MACS1206 field, Gaussian smoothed to the 
same resolution of 6*fwhm = 1-5'. The middle and right 
panels display the number and (i^-corrected) i? c -band 
luminosity density fields, respectively, of green cluster 
galaxies (see Table 3). 

Overall, mass and light are similarly distributed in the 
cluster: The cluster is fairly centrally concentrated in 
projection, and associated with elongated LSS running 
north-west south-east (NW-SE), both in the projected 
mass and galaxy distributions. A more quantitative char- 
acterization of the 2D matter distribution around the 
cluster will be given in Section 6. 

5.2. Cluster Weak-Lensing Profiles 



12 



CLASH: Full Lcnsing Analysis of MACS1206 



Weak lensing convergence K=£/£ cr 



E 
o 

s. 

o 

<D 

D 

a> 

o 





-o.\ 


, 


0.1 


0.2 1 


0.3^ 




5 -5-10 
Offset RA (arcmin) 



Galaxy number density, £„ (arcmin 
o 



Luminosity density, E, (h 10 L sun /arcmin ) 
o 




5 0-5 
Offset RA (arcmin) 



5 0-5 -10 
Offset RA (arcmin) 



Figure 8. Comparison of the surface-mass density field and the cluster galaxy distributions in MACS1206. Left: linear reconstruction 
of the dimcnsionlcss surface-mass density field, or the lensing convergence k(0) = E(0)/S cr ; t , reconstructed from Subaru distortion data. 
Middle: observed surface number density distribution S n (0) of green galaxies, representing cluster member galaxies. Right: observed R c - 
band surface luminosity density distribution S; (0) of the same cluster membership. The solid ellipse in each panel indicates the respective 
mean cllipticity and orientation measured within a circular aperture of 8', which is slightly larger than the cluster virial radius (6* v ; r ss 6.9'). 
The pair of gray solid lines in each panel defines the north-west (NW) and south-east (SE) excess regions. All images are smoothed with 
a circular Gaussian of FWHM 1.5'. The field size is 24' X 24'. North is to the top, east to the left. 



Now we derive azimuthally-averaged lens distortion 
and magnification profiles from the Subaru data. We cal- 
culate the weak-lensing profiles in N discrete radial bins 
from the cluster center (Section 4.2), spanning the range 
[$min , $max] with a constant logarithmic radial spacing 
Alnf9 = ln(6> max /f5 m j n )/Af, where the inner radial bound- 



ary 



is taken to be 9 



lin = 0.8' (> Bin ). The 
outer radial boundary # max is chosen to be 6> max = 16' 
(fl max rj 3.8 Mpc h^ 1 ), sufficiently larger than the typical 
virial radius r v ; r of high mass clusters (r v ; r « 1.6 Mpc /i -1 
for MACS1206; see Section 6), but sufficiently small with 
respect to the size of the Suprimc-Cam's ficld-of-vicw so 
as to ensure accurate PSF anisotropy correction. The 
number of radial bins is set to N = 8, chosen such that 
the detection S/N (defined as in Equation 12) is of the 
order of unity per pixel. 

5.2.1. Lens Distortion 

For each galaxy, we define the tangential distortion 
g+ and the 45° rotated component, with respect to the 
cluster center, from linear combinations of the distor- 
tion coefficients (c/i, (72) as </+ = —(171 cos 2<j) + g% sin 2(f)) 
and g x = — (<?2 cos 2r^ — g\ sin 20), with <fi being the po- 
sition angle of an object with respect to the cluster cen- 
ter. In the absence of higher-order effects, weak lensing 
only induces curl-free tangential distortions, while the 
azimuthal averaged x component is expected to vanish. 
In practice, the presence of x modes can be used to check 
for systematic errors. 

For each galaxy sample, we calculate the weighted av- 
erage of g+ in a set of radial bins (i = 1, 2, N) as 



9+(0i 
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(13) 



where the index k runs over all objects located within 
the ith annulus, 0* is the weighted center of the ith ra- 
dial bin, and the weight factor w^) is defined by Equa- 
tion 9. We use the continuous limit of the area-weighted 
center for 0i (see Appendix A of Umetsu & Broadhurst 
2008). We perform a bootstrap error analysis to assess 
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Figure 9. Cluster weak-lensing radial profiles as measured from 
background galaxies registered in deep Subaru images. The top 
panel shows the tangential reduced shear profile g+(9) (squares) 
based on Subaru distortion data of the full background (red+blue) 
sample. The bottom panel shows the count depletion profiles n(8) 
due to magnification for a flux-limited sample of red background 
galaxies. The circles and triangles show the respective results with 
and without the mask correction due to bright foreground objects 
and cluster members. The horizontal bar represents the constraints 
on the unlensed count normalization, no, as estimated from Subaru 
data. Also shown in each panel is the joint Bayesian fit (68% CL) 
to both profiles. 

the uncertainty er +! i in the tangential distortion profile 
g + s (Umetsu et al. 2010). 

In Figure 5, we compare azimuthally-averaged radial 
profiles of g + and g x as measured from our red, blue, 
green, and blue+red galaxy samples (Section 3.4). For 
all samples, the x component is consistent with a null 
detection well within 2cr at all radii, indicating the relia- 
bility of our distortion analysis. The red and blue popu- 
lations show a very similar form of the radial g + profile 
which declines smoothly from the cluster center. The 
observed tangential distortion signal is significant with a 
total detection S/N of 8.1 and 5.1 for the red and the 
blue sample, respectively, both remaining positive to the 
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limit of our data, (9 max = 16'. The detection significance 
is improved to 9.3a using a full composite sample of Sub- 
aru blue+red background galaxies (see the top panel of 
Figure 9). 

In Figure 5 we also compare the Subaru data with 
the results obtained from CFHT/Mcgacam data (Section 
3.1) using the same analysis pipeline as described in Sec- 
tion 3. For this we identified 15875 background galaxies 
(n g ~ 4.5 galaxies arcmin -2 ) with Megacam grz photom- 
etry using our CC background selection method (Section 

3.4) , and estimated a mean depth of z SiB s ~ 1.09, com- 
parable to that of the Subaru full background sample 
(zs ef£ = 1.15 ± 0.1, n g m 13 galaxies arcmin -2 ; Section 

3.5) . This comparison shows excellent agreement where 
the data overlap, demonstrating the robustness of our 
analysis. 

5.2.2. Magnification Bias 

For the number counts to measure magnification, wc 
follow the prescription of Umetsu et al. (2011b). We 
use a sample of red galaxies (Section 3.4), for which 
the intrinsic count slope s at faint magnitudes is rel- 
atively flat, s ~ 0.1, so that a net count depletion 
results (Broadhurst et al. 2005a; Umetsu & Broadhurst 
2008; Umetsu et al. 2010, 2011b). The blue background 
population, on the other hand, tends to have a steeper 
intrinsic count slope close to the lensing invariant slope 
(s = 0.4). 

The count-in-cell statistic N(9) is measured from a 
flux-limited sample of red background galaxies on a reg- 
ular grid of equal-area cells, each with a constant solid 
angle Af2. The practical difficulty here is contamination 
due to the intrinsic clustering of background galaxies, 
which locally can be larger than the lensing-induced sig- 
nal in a given cell. In order to obtain a clean measure of 
the lensing signal, such intrinsic clustering needs to be 
down- weighted and averaged over (e.g., Broadhurst et al. 
1995; Umetsu & Broadhurst 2008). 

To overcome this problem, we azimuthally average the 
red galaxy counts N(9) and obtain the average surface 
number density n„i = ria(9i) — (dN(6j,)/dQ) as a func- 
tion of radius from the cluster center (i = 1,2, N). 
Here we use the approach developed in Umetsu et al. 
(2011b) to account and correct for the masking effect 
due to bright cluster galaxies, foreground objects, and 
saturated objects. The errors a^.i for n^.i include both 
contributions from Poisson errors in the counts and 
contamination due to intrinsic clustering of red back- 
ground galaxies. Thanks to the wide field of view 
Subaru/Suprime-Cam, the normalization and slope pa- 
rameters for the red sample are reliably estimated as 
r;, = 11.4 ± 0.3 galaxies -2 and s = 0.133 ± 0.245 from 
the coverage-corrected source counts in the outer region 
( > 10'). 

We show in the bottom panel of Figure 9 the resulting 
magnification profile derived from our flux-limited sam- 
ple of red background galaxies (z' < 24.6 mag; see Table 
3). A strong depletion of the red galaxy counts is seen in 
the central, high-density region of the cluster and clearly 
detected out to ^ 4' from the cluster center. The statis- 
tical significance of the detection of the depletion signal 
is 4.4(7, which is about half the S/N ratio of the tangen- 
tial distortion derived from the full background sample 
shown in the top panel of Figure 9. The magnification 



measurements with and without the masking correction 
are roughly consistent with each other. 

5.3. Mass Profile Reconstruction 

The relation between observable distortion (<?) and un- 
derlying convergence (k) is nonlocal. Hence the mass dis- 
tribution derived from distortion data alone suffers from 
a mass-sheet degeneracy (§ 2). 

Here we construct a radial mass profile from com- 
plementary lens distortion and magnification measure- 
ments, {g+,i}f = i and following the Bayesian 
prescription given by Umetsu et al. (2011b), effectively 
breaking the mass-sheet degeneracy A brief summary of 
this Bayesian method is provided in Appendix A.l. The 
model is described by a vector s of parameters containing 
the discrete convergence profile {ni}f =1 in the subcritical 
regime (0, > #Ein)> and the average convergence within 
the inner radial boundary 9 m i n of the weak-lensing data, 
K m in = k(< 6* m in), so that s = {7t min , n t }f =1 , being spec- 
ified by (N + 1) parameters. 

Wc find a consistent mass profile solution s based on a 
joint Bayesian fit to the observed distortion and magni- 
fication measurements, as shown in Figure 9. The detec- 
tion significance has been improved from 9.3er to 11.4a by 
adding the magnification measurement, corresponding to 
an improvement by ~ 23%, compared to the lensing dis- 
tortion signal (Umetsu et al. 2011b; Coe et al. 2012). 

The resulting mass profile s is shown in Figure 6, along 
with our primary strong- lensing model (Sections 4.1- 
4.3). Our independent strong- and weak-lensing mass 
profiles are in good agreement where they overlap, and 
together form a well-defined mass profile. The outer mass 
profile derived from weak lensing exhibits a fairly shallow 
radial trend with a nearly isothermal logarithmic den- 
sity slope in projection, 72D = — rflnS/rflni? ~ 1. Note, 
this flat behavior is not clearly evident in the tangential 
distortion profile, which is insensitive to sheet-like mass 
overdensitics (Section 2). To constrain the cluster prop- 
erties from the composite halo+LSS mass profile, this 
LSS contribution needs to be taken into account and cor- 
rected for. We will come back to this point in Sections 
6.2 and 6.4. 

Also shown in Figures 6 and 7 is a purely shear-based 
reconstruction using the one-dimensional (ID) method 
of Umetsu & Broadhurst (2008, see also Umetsu et al. 
2010), based on the nonlinear extension of aperture mass 
densitometry (Clowe et al. 2000). Here we have adopted 
a zero-density boundary condition in the outermost ra- 
dial bin, 16' < 9 < 18'. The total S/N ratio in the recov- 
ered mass profile is w 9.2, which agrees well with rj 9.3 
in the g + profile (Section 5.2.1). Our results with differ- 
ent combinations of lensing measurements and boundary 
conditions, having different systematics, are in agreement 
with each other. This consistency demonstrates that our 
results are robust and insensitive to systematic errors. 

The projected cumulative mass profile M 2 d(< 0) is 
given by integrating the density profile s = {7« m i n , «n}i=i 
(see Appendices A and B of Umetsu et al. 2011b) as 

M 2D (< 6i) = 4A0) 2 £critK m i„+27rD 2 £ crit / d\n99 2 K (9). 

(14) 

We compare in Figure 7 the resulting M 2 d profiles dc- 



14 



CLASH: Full Lcnsing Analysis of MACS1206 



rived here from a wide variety of strong- (Section 4) and 
wcak-lcnsing analyses, along with the model-independent 
Einstein-radius constraints of M2d(< 17") = 5. x 
10 13 Af Q h- 1 at fen = 17" ± 2" (z s = 1.03) and M 2D (< 
28") = l.ltai x 10 14 M Q h- 1 at 6» Ein = 28" ± 3" (z s = 
2. 54). 44 Again, we find good agreement in the regions of 
overlap among the results obtained from a variety of lens- 
ing analyses, ensuring consistency of our lensing analysis 
and methods. 

Unlike the non-local distortion effect, the magnifica- 
tion falls off sharply with increasing distance from the 
cluster center. For MACS1206, we find k 1% at radii 
^,10', where the expected level of the depletion signal is 
n^/no — 1 ~ —2k for a maximally-depleted sample with 
s = 0, indicating a depletion signal of 2% in the outer 
region where we have estimated the unlensed background 
counts, no- This level of signal is smaller than the frac- 
tional uncertainties in estimated unlensed counts no of 
3% (Section 5.2.2), thus consistent with the assumption. 
Note that the calibration uncertainties in our observa- 
tional parameters (no, s, w) have been marginalized over 
in our Bayesian analysis (Appendix A). 

In the presence of magnification, one probes the num- 
ber counts at an effectively fainter limiting magnitude: 
nicut + 2.5 log 10 fi(9). The level of magnification is on 
average small in the wcak-lensing regime but reaches 
fi f* 1.6 (at z s e 1.1) for the innermost bin in this clus- 
ter. Hence, we have implicitly assumed in our analysis 
that the power-law behavior (Equation [7]) persists down 
to ~ 0.5 mag fainter than m cu t where the count slope 
may be shallower. For a given level of count depletion, 
an underestimation of the effective count slope could lead 
to an underestimation of /x, thus biasing the resulting 
mass profile. However, the count slope for our data flat- 
tens only slowly with depth varying from s ~ 0.13 to 
~ 0.05 from a limit of z' = 24.6 to 25.1 mag, so that 
this introduces a small correction of only ~ 10% for the 
most magnified bins (/i ~ 2). In fact, we have found a 
good consistency between the results with and without 
the magnification data. 

6. MASS PROFILE FROM JOINT WEAK AND STRONG 
LENSING ANALYSIS 

In this section, we aim to quantify and characterize the 
mass distribution of MACS1206 using our comprehen- 
sive lensing measurements derived from the deep HST 
and Subaru observations described in Sections 4 and 5. 
Here, we compare the cluster lensing profiles with the 
theoretically and observationally motivated NFW model 
(Navarro et al. 1997) to characterize the cluster mass 
profile. Our use of the NFW model enables the most 
direct comparison with detailed theoretical predictions 
for the internal structure of DM halos based on A-body 
simulations (e.g., Duffy et al. 2008; Klypin et al. 2011; 
Prada et al. 2011; Bhattacharya et al. 2011). The choice 
of profile shape does not significantly affect the derived 
halo concentrations (e.g., Duffy et al. 2008). 

To be able to constrain the inner density slope, we con- 
sider a generalized parameterization of the NFW model 

44 Zitrin et al. (2012) quote their full-model based estimates on 
the respective integrated masses of M20 = 6 ± 0.7 X 10 13 Mq h~ 1 
and A/ 2D = 0.94 ± 0.11 x 10 14 M o h,- 1 . 



(gNFW, hereafter) of the form (Zhao 1996; Jing & Suto 
2000): 

{r/r s ) a (l + r/r s y 

where p s is the characteristic density, r s is the character- 
istic scale radius, and a represents the inner slope of 
the density profile. This reduces to the NFW model 
for a = 1. We introduce the radius r_2 at which 
the logarithmic slope of the density is isothermal, i.e., 
73D = 2. For the gNFW profile, r_ 2 = (2 — a)r s , and 
thus the corresponding concentration parameter reduces 
to c_ 2 = r V i r /r_2 = c vir /(2 — a). We specify the gNFW 
model with the central cusp slope, a, the halo virial mass, 
Mvir, and the concentration, c_2 = c v i r /(2 — a). We em- 
ploy the radial dependence of the gNFW lensing profiles 
given by Keeton (2001). 

6.1. Model- Independent Constraints 

First, we constrain the NFW model parameters 
p = (Af v i r , c v i r ) by combining model-independent weak- 
lcnsing distortion, magnification, and strong-lensing 
Einstein-radius measurements, whose systematic errors 
are well understood from numerical simulations (e.g., 
Meneghetti et al. 2011; Rasia et al. 2012). The x 2 func- 
tion for the combined Einstein-radius and weak-lensing 
constraints is expressed as 

X 2 =x!in+XWL> ( 16 ) 

where the xlm f° r t nc Einstein-radius constraints is de- 
fined by (see Umetsu & Broadhurst 2008; Umetsu et al. 
2010) 

2 _ V- [1 -g+,i(P.Z»,«)] / 17 N 
AEin — 2-^i 2 \ Ll ) 

i=l (7 +- i 

with Agin being the number of independent Einstein- 
radius constraints {#Ein,;}iIfi n from sources with differ- 
ent redshifts {z s ,i}^i and g+,i(p, z s>i ) = g(9 E m,i\p, z s ,i) 
being the NFW model prediction for the reduced tangen- 
tial shear at 8 = 6*Em,i, evaluated at the source redshift 
z s = z s< i. Note, the Einstein radius marks the point of 
maximum distortion, g + = (k — k)/(1 — k) = 1: i.e., 
K = 1 within #Ein- The x 2 function for our full weak- 
lensing analysis (Section 5.3) is described by 

(18) 

where s = {n m i n , Ki}fL 1 is the mass profile reconstructed 
from the combined lens distortion and magnification 
measurements, s(p, z Sj0 ff) is the NFW model prediction 
for s, and C\vl is the full covariance matrix of s defined 

as 

C WL =C+C lss (19) 

with C being responsible for statistical measurement er- 
rors (Appendix A.l) and C lss being the cosmic covari- 
ance matrix responsible for the effect of uncorrelated LSS 
along the line of sight (Hockstra 2003; Hockstra et al. 
2011; Umetsu et al. 2011a; Oguri & Takada 2011). 45 In 

45 As discussed in Oguri et al. (2010),for a ground-based weak- 
lensing analysis, the shot noise is a more dominant source of the 
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all modeling below, the effective source redshift z s , e ff = 
1.15 ± 0.1 of our full background sample is treated as 
a nuisance parameter, and its uncertainty is marginal- 
ized over. In order to evaluate C lss , we assume the con- 
cordance ACDM cosmological model of Komatsu et al. 
(2011) and use the fitting formula of Peacock & Dodds 
(1996) to compute the nonlinear matter power spectrum. 
We project the matter spectrum out to an effective source 
redshift of z s>e s = 1-15 to calculate C lss for weak-lensing 
observations. For details, see Umetsu et al. (2011a). For 
Einstein-radius measurements, we conservatively assume 
an rms displacement of 2" due to uncorrelated LSS, as 
predicted by recent theoretical work (~ 2" for a distant 
source at z s ~ 2.5; see Host 2012; Julio et al. 2010). This 
is combined in quadrature with the measurement error 
in #Ein (Table 1) to estimate a total uncertainty 0+^. 46 

For strong lensing, we use double Einstein-radius con- 
straints (A^Ein = 2) from the multiple image systems at 
z s = 1.03 and z s = 2.54 (Table 1). For weak lensing, 
the cluster mass profile s is measured in N + 1 = 9 bins. 
Hence, we have a total of 11 constraints. 

The resulting constraints on the NFW model parame- 
ters are summarized in Table 5. 

6.1.1. Weak Lensing Constraints 

First of all, when no magnification or strong-lensing in- 
formation is included, the best-fit model is obtained from 
a tangential reduced shear fitting as M v ; r = 0.99^Q 2g x 

10 15 M Q h- 1 and c vir = 5.7±|:f with Xmin/dof = 3.3/6. 47 
Next, when magnification bias is included to break 
the mass-sheet degeneracy, we find M v ; r = LISIq ^ x 

IO^Mq/j" 1 and c vir _= 4.0±?;J (x^ in /dof =_ 4.5/7), 
which is consistent within the large uncertainties with 
the purely shear-based results, but is in favor of a larger 
M v j r and a smaller c v i r , owing to the shallow outer mass 
profile. This is demonstrated in the bottom panel of 
Figure 6, which shows significant deviations An from 
our reference NFW model (M vir 1.1 x 1O 15 M hr x 
and Cvir ~ 6.9; see Section 6.4) at cluster outskirts, 
R ^ 1 Mpc h^ 1 (9 £ 4'). This large-scale excess in pro- 
jected mass is also shown in Figure 7 in terms of the 
integrated projected mass profile M 2 d(< R)- Both fits 
here underestimate the observed Einstein radius (see Ta- 
ble 5). 

6.1.2. Combining Einstein Radius Constraints with Weak 

Lensing 

When the inner Einstein-radius information is com- 
bined with weak lensing, we obtain tighter parameter 
constraints. By combining all lens distortion, magnifica- 
tion, and Einstein-radius constraints (Equation 16), we 
find M vil = 1.0±g | x 10 15 M Q /i- 1 and c vir = 6.8^ 

measurement errors than the cosmic noise contamination. They 
found from a weak lensing analysis of 25 X-ray selected clusters 
that the best-fit parameters are not largely biased by including the 
cosmic noise covariance, but arc in general consistent with each 
other within statistical uncertainties. 

46 Following Umetsu & Broadhurst (2008), we propagate the 
uncertainty in $Ein to assuming a singular isothermal sphere 
(SIS) model. At r^r s , the density slope of NFW is shallower 
than that of SIS (see Figure 1 of Wright & Brainerd 2000), so that 
this gives a fairly conservative estimate of CT+(#Ein)- 

47 We follow Hockstra (2003) to calculate the cosmic shear co- 
variance matrix. 



(Xmin/dof = 6.9/9), corresponding to an effective Ein- 
stein radius of 6*Ei n ~ 26" at z s = 2.5. That is, a slightly 
higher concentration is favored to reproduce the observed 
large Einstein radii (Broadhurst & Barkana 2008). 

6.2. Mass and Galaxy Distribution Shapes in and 
around the Cluster 




Figure 10. Logarithmically-scaled XMM-Newton mosaic, 

exposure-corrected image of MACS1206 in the 0.5-2 keV band, 
smoothed with a Gaussian of 8.0" FWHM. Overlaid are contours 
from the exposure-corrected Chandra 0.5-2 keV image, smoothed 
with a Gaussian of 1.5" FWHM. The field size is 7.5' X 6.0', with 
north to the top and east to the left. The scale bar shows 2' or 
about 680 kpc « 1.1»*2500- X-ray emission is concentrated around 
and peaked on the BCG, but shows some elongation within 9 ;$ 1' 
at a position angle around 120 deg east of north, aligned with the 
orientation of the projected mass distribution. At larger distances 
from the cluster center, the cluster appears fairly round in both 
Chandra and XMM images. 

The presence of surrounding large scale structure in 
MACS1206 has a non-negligible impact on the deter- 
mination of cluster mass profile especially at large radii 
(Sections 5.3 and 6.1). It is therefore necessary to assess 
and correct for their effects on the projected mass profile. 
Here we use two different methods to quantify the ellip- 
ticity and orientation of the projected mass distribution 
in and around the cluster. 

First, following the prescription given by Umetsu et al. 
(2009), we introduce mass-weighted quadrupole shape 
moments around the cluster center, in analogy to Equa- 
tion (8), defined as 

Q a p= [ d 2 6A6 a A6pZ(0) (a,/3 = 1,2),(20) 

JA6»<e max 

where 9 max is the circular aperture radius, and A9 a is 
the angular displacement vector from the cluster center. 
We construct with {Q a p} a spin-2 ellipticity measure 
es = \e^\e 2l ^ c , where the ellipticity is defined such that, 
for an ellipse with major and minor axes a and b, it re- 
duces to | eg | = 1 — b/a and <p e is the position angle of the 
major axis (Bertin & Arnouts 1996), measured north of 
west here. Similarly, the spin-2 ellipticity for the cluster 
galaxies is defined using the surface number and i? c -band 
luminosity density fields of CC-selectcd cluster galaxies 
(Section 5.1), £„(0) and £;(#). We calculate weighted 
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moments using only those pixels above the 2a threshold 
with respect to the background level (estimated with the 
biweight scale and location, see Beers et al. 1990). Prac- 
tical shape measurements are done using pixelized maps 
shown in Figure 8. 

Next, we constrain the ellipticity and orientation of the 
projected mass distribution by directly fitting a 2D shear 
map with a single elliptical lens model. Here, we closely 
follow the prescription given by Oguri et al. (2010) to 
construct an elliptical NFW (eNFW, hereafter) model 
(see also Oguri et al. 2012), by introducing the mass el- 
lipticity |e E | = 1 — b/a in the isodensity contours of the 
projected NFW profile £(i?) as R 2 -> X 2 (l - |e E |) + 
Y 2 '/{1 - |e s |) (for details, see Oguri et al. 20 1 0). 48 The 
model shear field is computed by solving the 2D Poisson 
equation (Keeton 2001). We then construct from Subaru 
data a lens distortion map (gi(O),g2(0)) and its covari- 
ance matrix C g (Equation A5) on a 2D Cartesian grid 
with V spacing, centered at the BCG. We exclude from 
our analysis the five innermost cells lying in the central 
region, 9 < 1', to avoid systematic errors (see Appendix 
A. 2). The halo centroid is fixed to the BCG position. 
Accordingly, the eNFW model is specified by four model 
parameters, p = (M V i r ,c V ir, |es|,</> e ). The constraints on 
individual parameters are obtained by projecting the 2D 
shear likelihood function (Equation (A6) in Appendix 
A. 2) to the parameter space (or, minimizing \ 2 )- 

In Table 6, we summarize our cluster ellipticity and ori- 
entation measurements. In this analysis, we are mainly 
interested in the orientation of the ellipticity, in order to 
correct for the effects of LSS along the axis of elonga- 
tion. An overall agreement is found between the shapes 
of mass, light, and galaxy distributions in MACS1206, 
especially in terms of orientation (Figure 8), within 
large uncertainties (Table 6). The mass distribution in 
and around the cluster is aligned well with the lumi- 
nous galaxies in the green sample, comprising mostly 
of cluster member galaxies (Section 3.4). For all cases, 
the position angle <j> e of the major axis is found to be 
fairly constant with radius # m ax, and lies in the range 
15 deg < <t> e < 30 deg. 

In the central region, we find a projected mass elliptic- 
ity of |es| ~ 0.3 and a position angle of e ~ 14 deg from 
the Pixelens analysis; we obtain consistent values for 
both |es| and cj> e from a different strong- lensing analysis 
(C. Grillo et al., in preparation) using only System 7 of 
Zitrin et al. (2012). A similar value is found for the pro- 
jected mass ellipticity of \e^\ — 0.26 ±0.16 (4> e ~ 19 deg) 
at #max = 4' using the weak-lcnsing £ map. From an 
elliptical King model fit to Chandra X-ray data (Figure 
10; for details, see Section 7.4), we find an ellipticity 
of 0.30 ± 0.03 (a/b « 1.5) and <f> e = 21.9 ± 1.7 deg at 

$max = 1.5 . 

On the other hand, we obtained higher values of ellip- 
ticity on large angular scales beyond the cluster virial 
radius, 6» vir = r vir /A ~ 7'. We find |e E | - 0.4-0.5 
at # m ax = 8' using the pixelized cluster mass, galaxy, 
and light distributions. From the 2D shear fitting to 
a single eNFW model, the projected mass ellipticity is 

48 As noted by Oguri et al. (2010), this elliptical model includes 
a triaxial halo model (Jing & Suto 2000) which gives a better de- 
scription of CDM halos in TV-body simulations than the spherical 
model. 



constrained in the range |es| = 0.68_o' 2 g (l e s| ~ 0.4 or 
a/b}^ 1.7 at la) at 6> max = 8'. This apparent increase 
in ellipticity with radius could be partly explained by 
the additional contribution from the surrounding LSS 
that is extended along the cluster major axis. Note, the 
observed tendency for the shear-based method to yield 
higher ellipticity estimates, compared to the mass-map 
based method, could be due to the nonlocal nature of the 
shear field, in conjunction with our single-component as- 
sumption in the 2D shear fitting analysis. Overall, this 
level of ellipticity is consistent within large errors with 
the mean cluster ellipticity (|es|) = 0.46 ±0.04 obtained 
by Oguri et al. (2010) from a 2D weak-lensing analysis 
of 25 X-ray luminous clusters. 

In what follows, we fix the position angle of the NW- 
SE cluster-LSS major axis to a reference value of 4> e = 
20 deg, which is close to the values derived from the 
Chandra X-ray data, S; and £ maps. We note that, 
in principle, the X-ray structure in a triaxial system 
is expected to be tilted with respect to the total mat- 
ter in projection, even in the absence of intrinsic mis- 
alignments (see Romanowsky & Kochanek 1998). In the 
present work, we define the NW and SE excess regions re- 
spectively as NW and SE outer cone regions with 9 > 4' 
centered on the cluster center, with opening angle 90 deg 
and position angle <f> e = 20 deg, as defined by the pair of 
gray solid lines in each panel of Figure 8. 

6.3. BCG-Cluster Alignment 















• 
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Figure 11. Detailed model fits to the BCG and its nine nearby 
galaxies in the ACS F814 image (« 50" X 45"). The panels show 
the ACS image (left), best-fit model (middle), and image residuals 
(right) after subtraction of the model. No systematic deviations are 
seen in the residuals between the data and the model, suggesting 
the BCG has not undergone any major merger recently. North to 
the top, east to the left. 

We have also obtained CLASH constraints on the mean 
BCG ellipticity and position angle derived from the ACS 
F814W image. For this we performed a detailed struc- 
tural analysis on the BCG using the snuc task in the 
X Vista software package. 49 In Figure 11 we show the 
ACS F814W image, best-fit model, and image residu- 
als after subtraction of the model. No systematic de- 
viations are seen in the residuals between the data and 
the model, suggesting the BCG has not undergone any 
major merger recently. The radial profiles of elliptic- 
ity and position angle were measured in several inde- 
pendent radial bins (0.2" ^ 9 ^ 10"), and their respec- 
tive (sensitivity weighted) mean values were obtained as 
(e s ) = 0.53 ± 0.03 and 4> e = (15.0 ± 2.3) deg (Table 6). 
Consistent results were found in several other HST bands 
(ACS F475W to F814W and WFC3 F105W to F160W). 
The mean BCG ellipticity is found to lie in the range 

49 http:/ /astronomy.nmsu.edu/holtz/xvista/index.html 
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0.46-0.53 with a small scatter of 0.02 across the ACS and 
WFC3 bands. The BCG position angle is constrained to 
be 4> e = (15.2 ± 0.4) deg, which is in excellent agreement 
especially with that derived independently from the large 
scale distribution S n of galaxies. 

6.4. Effects of Surrounding Large Scale Structure 

In this subsection, we look into the azimuthal depen- 
dence of the radial projected mass distribution, </>), 
to assess and correct for the effect of surrounding LSS 
on the cluster mass profile measurement. Because of the 
non-local nature and inherent insensitivity to sheet-like 
overdensitics of the shear field, it is essential to use the 
combination of lens magnification and distortion to re- 
construct the projected cluster mass distribution embed- 
ded in large scale structure. For this purpose, we extend 
the ID Bayesian method of Umetsu et al. (2011b) into a 
2D mass distribution by combining the 2D shear pattern 
g(6) with the azimuthally-averaged magnification mea- 
surements n M (0). In the 2D analysis, our model s is a 
vector of parameters containing a set of discrete mass ele- 
ments on a grid of N cc \\ independent cells, s = {/«m}m=i- 
A brief summary of this 2D method is given in Appendix 
A. 2. The details of the method will be presented in our 
forthcoming paper (K. Umetsu et al., in preparation). 

By combining Subaru distortion and magnification 
data, we construct here a mass map over a 30 x 30 grid 
with 0.8' spacing, covering a 24' x 24' field around the 
cluster (iVceij = 900). We have 2 x 896 distortion con- 



straints {<?i(0m)} m =i and {g2{0 m )} r n=[ over the mass 
grid, excluding the four innermost cells lying in the clus- 
ter central region (see Appendix A. 2), and N = 8 radial 
magnification constraints {n^i{9i)}f =1 . Hence, we have a 
total of 1800 constraints (900 degrees of freedom). Addi- 
tionally, we marginalize over the calibration uncertainties 
in the observational parameters (no,s,u>) (Section 3.5). 
The best solution s has been obtained with x^ in /dof = 
1058/900. We then follow Umetsu & Broadhurst (2008) 
to calculate the radial mass distribution (S(i?)) and its 
covariance matrix from a weighted projection of the k 
map, where we conservatively limit our 2D analysis to 
radii smaller than = 12' (R « 2.9 Mpc /i" 1 ). We check 
that the azimuthally-averaged radial mass profile con- 
structed from the k map reproduces our corresponding 
ID results (Section 5.3). 

We show in Figure 6 the radial mass distribution ob- 
tained excluding the NW and SE excess regions (defined 
in Section 6.4; see also Figure 8). This weak-lensing mass 
profile, corrected for the effect of surrounding LSS, ex- 
hibits a steeper radial trend than that averaged over all 
azimuthal angles. We note, a slight remaining excess 
is seen at 8 ^ 5' (R £ 1.2 Mpc ft -1 ), which may be as- 
sociated with the likely north-south filamentary feature 
(Section A. 2). By fitting the "LSS-corrected" mass pro- 
file with an NFW profile, we find a higher concentra- 
tion C v ir 

,2 



7.5±?;| with M v j r = 1.15;^2o x 1O 15 M h~ 



(Xmin/dof = 10.6/6). This model predicts an Einstein 
radius of #Ein ~ 32" for z s = 2.5, comparable to the 
observed value, Ein = 28" ± 3". 

6.5. Full Lensing Constraints 

As shown in Figures 6 and 7, our weak and strong lens- 
ing data agree well in their region of overlap. Here we 
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Figure 12. Constraints on the NFW model parameters 
(M v i r , c v i r ), the halo virial mass and concentration, derived from 
weak-lensing (red), strong- lensing (blue), and joint weak and strong 
lensing (yellow) analyses. The weak-lensing results are obtained 
with the LSS correction (Section 6.4). The contours show the 
68.3% and 95.4% confidence levels, estimated from Ax 2 = x 2 ~ 
xtiir, = 2.3 and 6.17, respectively. The circles indicate the re- 
spective best-fit model parameters. For weak lensing, the source 
redshift uncertainty, z s = 1.15 ± 0.1, is marginalized over. 
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Figure 13. Constraints on the gNFW model parameters, namely, 
the central cusp slope a, the halo virial mass M v - lr , and the halo 
concentration c_2 = r v j r /r_2 = c v ; r /(2 — a), when all of them 
are allowed to vary, derived from combined weak and strong lens- 
ing. The weak-lensing results arc obtained with the LSS correction 
(Section 6.4). The left and right panels show the two-dimensional 
marginalized constraints on (M v ; r , a) and (c_2 , a), respectively. In 
each panel of the figure, the contours show the 68.3%, 95.4%, and 
99.7% confidence levels, and the circle indicates the best-fit model 
parameters. 

further improve the statistical constraints on the halo 
parameters p = (M v - ir , c v i r , ol) by combining the joint 
weak-lensing distortion and magnification constraints 
Xwl(p> z s,ctf) (Section 6.1) with the inner mass pro- 
file Ki based on the detailed strong-lensing analysis of 
Zitrin et al. (2012). 

We write the combined x 2 function of our full lensing 
constraints as 



X 



XWL 



XSL 



with Xsl f° r strong lensing being defined as 



Xsl 



J2 [*i - iP)} (Csl)" 1 [ Kj - ^ (p)} , (22) 



where is defined in 25 discrete bins over the radial 
range [4", 53"] (see Section 4) and scaled to a fidu- 
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cial depth z s = 2.54 of the strong-lensing observations, 
matched to the spectroscopically-confirmed five image 
system (System 4 of Zitrin et al. 2012), hi is the theo- 
retical prediction for m\ Csl = C + C\ ss is the bin-to-bin 
covariance matrix for the discrete k profile, with C being 
the self-calibrated covariance matrix derived in Section 
4.3 and Ci ss being the cosmic noise contribution. We use 
a consistent single source plane at z s = 2.54 to evaluate 
C lss . 

The resulting NFW and gNFW fits are summarized 
in Tables 5 and 7, respectively. For both models, 
we show the respective fits derived with and with- 
out the LSS correction for the outer weak-lcnsing pro- 
file (R > lMpc/i" 1 ). We find that, when the detailed 
strong-lensing information is combined with weak lens- 
ing, the LSS correction does not significantly affect 
the fitting results with the adopted NFW/gNFW form. 
Moreover, all these models properly reproduce the ob- 
served location of the Einstein radius, #Ein ~ 28". 

Here we summarize our primary results obtained with 
the LSS correction. The confidence contours on the 
NFW parameters (M v ; r ,c v i r ) are shown in Figure 12. 
The constraints are strongly degenerate when only the 
inner or outer mass profile is included in this fit. Com- 
bining complementary weak and strong lcnsing informa- 
tion significantly narrows down the statistical uncertain- 
ties on the NFW model parameters, placing tighter con- 
straints on the entire mass profile (Model #7 of Table 
5): M vir = imtoie x 10 15 Af o /i" 1 and c vir = 6.9±H 
with Xmin/dof = 18.0/31, corresponding to a Q value 
goodncss-of-fit of Q = 0.970. Next, when a is allowed to 
vary (Table 7), we find M vir = 1.06±g;f| x 10 15 M Q ft" 1 , 
c_2 = 7.0±\%, and a = 0.97±£g with Xmin/ dof = 
18.0/30 and Q = 0.960, being consistent with the sim- 
ple NFW form with a = 1. Thus, the addition of the 
a parameter has little effect on the fit, as shown by the 
quoted x 2 an d Q values. The two-dimensional marginal- 
ized constraints on (M v ; r ,a) and (c_2, Q() are shown in 
Figure 13. 

6.6. Impact of the Choice of Strong Lensing Models 

In this subsection we address the impact of the choice 
of strong-lensing models on the determination of the halo 
mass and concentration parameters in a joint weak and 
strong lensing analysis. As an alternative choice to the 
Zitrin et al. (2012) model, we consider here Pixelens 
(non-parametric) and Lenstool (parametric) models, 
in combinations with our LSS-corrected weak-lensing 
mass model (Section 6.4). For each case, we define the 
X 2 function for strong lensing as in Equation (22), and 
minimize the total x 2 function (Equation 21) to estimate 
the NFW parameters (M v i r ,c v i r ). 

The resulting model constraints are tabulated in Ta- 
ble 8. We find that both parameters based on differ- 
ent strong-lensing profiles are consistent with each other 
within the statistical errors. This also indicates consis- 
tency between these strong-lensing models and our weak- 
lcnsing measurements, as shown in Figure 7. We find a 
tendency for Pixelens to yield somewhat higher mass 
estimates compared to other strong-lens modeling meth- 
ods, as discussed by Grillo et al. (2010, their Appendix). 

When the NFW (gNFW) form is assumed, the 
Zitrin et al. (2012) model predicts a somewhat higher 



concentration and a lower mass than other models as 
implied by its correspondingly higher central density at 
^ 0.2' (see Figure 7). When the inner fitting radius 
is increased from 4" to 12" (~ 50kpc/i _1 ), we find a 
fractional increase of ~ 9% in M v ii and a fractional de- 
crease of - 17% in c vir (6.91q g -> 5. 7±\i). Including 
these variations as systematic uncertainties in our mass- 
concentration determination, the spherical NFW model 
for MACS1206 is constrained as M vir = (1.07±g;?g ± 
0.10) x IO^Mq/i- 1 and c vir = 6.9±J;!j ± 1.2, (statisti- 
cal followed by systematic uncertainty). Similarly, when 
the central 50kpc/i _1 region is excluded from the fit, 
we have Af vir = (1.17±g"JB ± 0.10) x 10 15 M Q h' 1 and 
Cvir = 5.T±\i ± 1.2. 

6.7. Alternative Mass Profile Fits 

Motivated by the apparently shallow projected density 
profile in the outer regions (cf. XMMU J2235. 3-2557 at 
z = 1.4, Jee et al. 2009), we consider here a softened 
power-law sphere (SPLS) model (Grogin & Narayan 
1996) as an alternative to the NFW profile, and perform 
profile fitting analyses on our full-range mass profile data 
(derived from Methods #6 and #7 in Table 5; see Sec- 
tions 6.4 and 6.5). 

The SPLS model has the same number of free param- 
eters as gNFW, namely three. The SPLS density profile 
is given by p(r) = po{l + r 2 /r 2 )(''~ 3 )/ 2 where po = p(0) is 
the central density, r c is the core radius, and the power- 
law index rj is restricted to lie in the range < rj < 2 
(Grogin & Narayan 1996). At r > r c , M(< r) cx r v . 
This reduces to a non-singular isothermal sphere (NIS) 
model when n — 1. The fitting results with and without 
the outer LSS correction (Methods #7 and #6, respec- 
tively) are summarized in Table 9. 

First, when rj is fixed to unity (NIS), the NIS model 
provides acceptable fits, but with larger residuals (x ) 
compared to the corresponding NFW fits with the same 
degrees of freedom (31): A X 2 = xLn,ms - xLlNFW = 
2.3 (Method #6) and 5.9 (Method #7) between the best- 
fit NIS and NFW models. Note, because of the asymp- 
totic M(< r) cx r behavior, the assumed NIS form leads 
to substantially higher masses at large radius (r 3> r c ) 
than what the NFW model predicts (~ 35% higher than 
the NFW values at r = 1.6Mpc/i _1 ). 

Next, when the outer slope is allowed to vary, the fit 
is noticeably improved for the results with the outer LSS 
correction (Method #7), corresponding to a difference 
of Ax 2 = xii„,Nis " xii„,sPLS = 4 - 4 between NIS and 
SPLS for 1 additional degree of freedom. For this, the 
best-fitting slope parameter is obtained as rj = 0.77lg 17 

(Xmin = 19-5 f° r 30dof), corresponding to 2.1 < 73d(^ 3> 
r c ) < 2.4. This SPLS model yields a virial mass of M V - 1T = 
(1.26 ±0.37) x lO 15 M /i- x (r vir w 1.73 Mpc^ 1 ). 

7. DISCUSSION 

7.1. Lensing Systematics 

Gravitational lensing probes the total mass projected 
on to the sky along the line of sight, so that the 
lensing-based cluster mass measurements are sensitive 
to projection effects arising from (1) additional mass 
overdensities (underdensitics) along the line of sight 
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(Meneghetti et al. 2010b; Rasia ct al. 2012) and (2) 
halo triaxiality (Hcnnawi et al. 2007; Oguri & Blandford 
2009; Meneghetti et al. 2010b; Rasia et al. 2012). 

7.1.1. Projection of Additional Mass Structures 

The first type of projection effects includes the cosmic 
noise from distant uncorrelated LSS projected along the 
same line of sight (Hockstra 2003) and massive structures 
within/around the cluster (i.e., cluster substructures and 
surrounding large scale filamentary structure). 

The former can not only increase statistical uncertain- 
ties but also produce covariance between radial bins. 
Accordingly this could bias the estimates of cluster pa- 
rameters. Our methods take into account the estimated 
contribution of cosmic covariance C\ ss in both weak and 
strong lensing profiles, and allow us to properly weight 
the weak and strong lensing when performing a combined 
halo fit. In our analysis, we find that the contribution of 
Ciss to the measurement errors is subdominant in both 
regimes; when the weak and strong lensing constraints 
are combined, the amount of degradation due to C\ ss is 
about 12% in the total S/N ratio. Thus the best-fit pa- 
rameters are not largely affected by including C\ ss , being 
consistent with each other within statistical uncertain- 
ties. 

The latter represents projection effects arising from 
the rich, substructured cluster environment. Recently 
Meneghetti ct al. (2010b) and Rasia ct al. (2012) used 
mock observations of simulated clusters in the ACDM 
cosmology to study the systematic effects in lensing and 
X-ray based mass measurements, finding that the stan- 
dard tangential-shear fitting method, assuming a sin- 
gle spherical NFW profile, can underestimate the true 
cluster mass Ma in the presence of massive substruc- 
tures, especially for low-mass systems. This is under- 
stood by noting that the azimuthally-averaged tangen- 
tial shear probes the differential surface mass density, 
7+(0) oc £(< R) - £(i?) (see Equation (5)). Rasia et al. 
(2012) found from their three most massive systems 
with M 2 oo > 7.5 x 1O 14 M h~ x that the level of bias 
is ~ —5% with no noticeable radial dependence at r = 
(^25007 r 1000, r 50o)- Our cluster mass estimate from the 
tangential-shear fitting is M vir = 0.99 to.26 x 10 15 M q h' 1 
(Model #1 of Table 5), which is about 7% lower than 
that from our NFW model based on the full-lcnsing con- 
straints (Model #7 of Table 5) from our comprehen- 
sive weak-lcnsing distortion, magnification, and strong- 
lcnsing analysis. This level of underestimation seems to 
be consistent with the simulation results of Rasia et al. 
(2012). 

7.1.2. Halo Triaxiality 

A degree of triaxiality is inevitable for collision- 
less gravitationally collapsed structures (Jing & Suto 
2000; Lemze et al. 2011), and can affect our cluster 
mass estimation (Oguri et al. 2005; Morandi et al. 2011; 
Sereno & Umetsu 2011; Sereno & Zitrin 2012). In the 
context of ACDM, prolate halo shapes are expected 
to develop by mass accretion along filaments at early 
stages of halo assembly; hence, dynamically-young, 
cluster-sized halos tend to have a prolate morphology 
(Shaw ct al. 2006; Lau et al. 2011). Accordingly, a large 
fraction of cluster-sized prolate halos, in the absence of 



selection bias, is expected to be elongated in the plane 
of the sky (Rasia ct al. 2012). On average, this will lead 
to an underestimation of the cluster mass in a statistical 
sense when a spherical deprojection (or forward mod- 
eling assuming a spherical halo) is applied (Rasia et al. 
2012). On the other hand, in the ACDM context, those 
clusters selected by the presence of giant arcs are likely 
to have their major axes closely aligned with the line of 
sight (Hennawi et al. 2007; Meneghetti et al. 2010b), be- 
cause this orientation boosts the projected surface mass 
density and hence the lensing signal. 

MACS1206 is an X-ray selected CLASH cluster 
(Postman et al. 2012), discovered in the MACS survey 
(Ebeling et al. 2001, 2009). For MACS1206, we find a 
large projected mass ellipticity of \e^\ = 1 — b/a £ 0.4 
(or a/b £ 1.7 at la) at large cluster radius (R £ r v ir « 
l^Mpc/i' 1 ) based on the Subaru weak-lensing analysis, 
where its position angle is well aligned with the BCG, op- 
tical, X-ray, and LSS shapes in projection space (Section 
6.4 and Table 6). The highly elliptical mass distribution 
in projection would suggest that its major axis is not 
far from the sky plane, and that its true mass and con- 
centration could be even higher than the projected mea- 
surements if the cluster size along the sight line is shorter 
than its effective size-scale (Vab) in the sky plane. 

7.2. Chandra and XMM-Newton X-ray Observations 
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Figure 14. Integrated total mass profiles M(< r) as a func- 
tion of spherical radius r derived from various observational probes 
(top). The blue shaded area shows the best-fit NFW model with 
la uncertainty from the combined weak and strong lensing mea- 
surements (Figures 6 and 7). The red solid lines represent the 
X-ray based NFW model (la confidence interval of the fit) derived 
using the JACO software from a simultaneous fit to Chandra and 
XMM- Newton observations. The enclosed masses based on Chan- 
dra data alone (solid line with error bars, orange) are derived as 
described in the text, assuming the parameterized pressure pro- 
file shape from Arnaud et al. (2010). The green square marks the 
Bolocam SZE mass estimate at the lensing-dcrived overdensity ra- 
dius r2500- The bottom panel shows the X-ray to lensing mass 
ratio Mx{< r)/M\ ena (< r) with la uncertainty as a function of 
radius r. The results are shown for both the Chandra-only and 
joint Chandra+XMM fits. 

Complementary multiwavclcngth observations serve as 
a useful guide to the likely degree of lensing bias. Here 
we retrieved and analyzed archival Chandra and XMM- 
Newton data of MACS 1206 to obtain an independent 
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Figure 15. Gas mass fraction profiles / gas (< r) = M gas (< 
r)/M(< r) as a function of spherical radius r derived from joint 
Subaru weak-lensing, Hubble strong-lensing, and Chandra X-ray 
observations. In each case the gas mass profile M gas (< r) is 
based on the Chandra X-ray data provided in the ACCEPT 
(Cavagnolo et al. 2009). The squares with error bars represent the 
results (Af g as/Mi eIls ) from the combined X-ray and lensing data 
without employing the hydrostatic equilibrium assumption. The 
circles with error bars show the Chandra-only results (M gas /Mx) 
based on the hydrostatic equilibrium assumption. The horizontal 
bar shows the constraints (68% CL) on the cosmic baryon fraction 
from the WMAP seven-year data, f b = n b /n m = 0.1675 ± 0.006. 

cluster mass estimate, as well as to constrain the physical 
properties of the X-ray gas. 

We perform a simultaneous fit to Chandra and XMM 
data-sets under the assumption that the intracluster gas 
is in hydrostatic equilibrium (HSE) with the overall clus- 
ter potential of the NFW form. The tool used for this 
analysis is JACO (Joint Analysis of Cluster Observa- 
tions, Mahdavi et al. 2007b); we refer the reader to this 
paper for the details of the X-ray analysis procedure, 
which we briefly summarize below. 

We use Chandra ObsID 3277 and XMM-Newton ob- 
servation 0502430401. We screen periods of flaring back- 
ground according to standard procedure, resulting in us- 
able exposure times of 23 ks and 26 ks, respectively. Ap- 
propriate coaddcd blank-sky fields allow us to subtract 
particle background spectra for both telescopes, and the 
residual (positive or negative) astrophysical background 
is included and marginalized over in the global cluster gas 
model. Spectra are extracted over seven annular bins for 
both Chandra and XMM-Newton. The extracted spec- 
tra extended out to a distance of 3.7' (1.26 Mpc), and 
contain an average of 1500 counts each. 

The model for the gas density distribution is a single 
/3-modcl multiplied by a power law of slope 7: 



Pair) = Po 



r 2 \ "30/2 

1 + 3 



(23) 



The power-law component is required to capture the 
steep increase of the density towards the center of the 
cluster; all parameters of the gas distribution are fit to 
the data. The metallicity is allowed to vary with ra- 
dius as well, as are the parameters of the NFW mass 
profile. Model spectra are generated self-consistently in 
concentric spherical shells and forward-projected onto 
the annular sky regions matching the extracted annuli. 



The resulting spectra are mixed using in-orbit energy- 
and position-dependent point spread functions for both 
Chandra and XMM-Newton. Systematic calibration un- 
certainties between Chandra and XMM-Newton spectra 
are taken into account by adding a 4% error (a typi- 
cal correction used in Mahdavi et al. 2008) in quadra- 
ture to each spectral bin used for the joint fits. This 
brings the joint % 2 into the acceptable range (x 2 = 1603 
for 1541 dof). An MCMC procedure is used to esti- 
mate errors on the best-fit quantities. After marginal- 
izing over all other parameters, we measure a total mass 
M 25 oo = (4.45 ± 0.28) x 10 14 M o , a gas mass M gas , 250 o = 
(0.54±0.02) x 10 14 M©, an NFW concentration parame- 
ter of C200 = 3.5 ± 0.5, an inner gas density profile slope 
of 0.7 ± 0.03, and a central cooling time of 2.1 ± 0.1 Gyr. 
In what follows, the examination of the X-ray results is 
conservatively limited to r < 1 Mpc. 

In Figure 14 wc plot the resulting X-ray based to- 
tal mass profile, M(< r), shown along with our NFW 
model from the full lensing analysis. The results of the 
NFW fit are also reported in Table 10. This X-ray model 
yields a total mass of M x = (4.6 ± 0.2) x 10 14 M Q at the 
lensing derived overdensity radius of r 2 5oo ~ 0.60 Mpc. 
This is in excellent agreement with the lensing mass 
at the same radius, Mi ons = (4.9 ± 0.9) x 10 14 M Q , 
which corresponds to the X-ray to lensing mass ratio, 
M x {< r 25 oo)/M lcns (< r 250 o) - ftn ^ +0 ' 23 



O2500 



0.95^. The 



(22500 value obtained here is in good agreement with re- 
sults from mock observations of 20 ACDM clusters by 
Rasia et al. (2012): a 250 o = 0.94±0.02. At this overden- 
sity, no significant bias was observed in detailed observa- 
tional studies by Zhang et al. (2008) and Mahdavi et al. 
(2008), who performed a systematic comparison of wcak- 
lcnsing and X-ray mass measurements for sizable cluster 
samples. In the bottom panel of Figure 14, we show the 
X-ray to lensing mass ratio oa as a function of cluster 
radius, in the radial range where X-ray observations are 
sufficiently sensitive. Overall, the mass ratio is consistent 
with unity especially at r ~ r 2 5oo- 

Ebeling et al. (2009) obtained a hydrostatic mass es- 
timate of M x = (1.7 ± 0.1) x 1O 15 M at r = 2.3 Mpc 
(their estimate for r 2 oo) assuming an isothermal (3 model 
with /3 = 0.57 ± 0.02 and their estimated temperature 
k B T = 11.6 ± 0.7keV in the radial range [70, 1000] kpc 
(Mx oc jS 1 ' 2 ! 1 ), which is high but consistent within the 
errors with M lons (< 2.3 Mpc) = (1.4 ± 0.3) x 1O 15 M 
obtained with our best NFW model based on the full 
lensing analysis. 

Our full lensing results, when combined with X-ray gas 
mass measurements (M gas ), yield a direct estimate for 
the cumulative gas mass fraction, / gas (< r) = M gas (< 
r)/M(< r), free from the HSE assumption. For this 
we use reduced Chandra X-ray data presented in the 
Archive of Chandra Cluster Entropy Profile Tables AC- 
CEPT (Cavagnolo et al. 2009). In Figure 15, we plot 
our / gas measurements as a function of cluster radius. 

We find a gas mass fraction of / gas (< r) = 13. 7j^Q% 
at a radius of r = 1 Mpc w 1.7r 2 5oo(~ 0.8rsoo), a typi- 
cal value observed for high mass clusters (Umetsu et al. 
2009; Zhang et al. 2009). When compared to the cos- 
mic baryon fraction /& = f2(,/Q m = 0.1675 ± 0.006 con- 
strained from the WMAP seven-year data (Jarosik et al. 
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2011), this indicates f gas /fb — 0.82+ ' 18 at this radius. 
At the innermost measurement radius r sa 40 kpc where 
the lensing and X-ray data overlap, we have / gas (< r) = 
3.4lJ'g%. Thus, the hot gas represent only a minor frac- 
tion of the total lensing mass near the cluster center, 
as found for other high-mass clusters (Lemze et al. 2008; 
Umetsu et al. 2009). 

Additionally, we derive a mass profile from sim- 
ulated annealing fits of the ACCEPT pressure pro- 
file (Cavagnolo et al. 2009), adopting the Arnaud et al. 
(2010, A10) "universal profile" (M. Donahue et al, in 
preparation). This Chandra-only mass profile is shown 
to be in good agreement with the lensing as well as joint 
Chandra+XMM results (Figure 14). The Chandra-only 
gas mass fraction profile is also shown in Figure 15. 

We conclude, on the basis of these results and com- 
parison with detailed statistical studies, that the level of 
orientation bias in this cluster is not significant given the 
large uncertainties in our lensing/X-ray observations, as 
well as the possible contribution from non-thermal pres- 
sure in the cluster core (e.g., Kawaharada et al. 2010). 

7.3. Bolocam SZE Observations 
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Figure 16. Bolocam SZE decrement images each 10' X 10', beam 
smoothed to an effective resolution of 82". Left: the processed im- 
age obtained when the data are filtered to remove atmospheric 
noise. The solid white contours denote S/N = —2, —4, —6, and 
the dashed contours denote S/N = +2, +4, .... Right: image ob- 
tained when the effects of the atmospheric filtering have been de- 
convolved to obtain an unbiased image of the cluster. 

We have also compared our lensing-derived results to 
mass estimates determined from the Sunyaev-Zcl'dovich 
effect (SZE) data. Using Bolocam at the Caltech Submil- 
limeter Observatory, we observed MACS1206 for approx- 
imately 11 hours in April 2011. These data were collected 
with Bolocam configured at an SZE-emission-weightcd 
band center of 140 GHz. Further details of the Bolocam 
instrument are given in Haig et al. (2004). We detect the 
cluster with an S/N value of 21.1 and a white noise rms 
of 24.9 /uKcMB-arcmin. We reduced these data accord- 
ing to the procedure described in detail in Sayers et al. 
(2011), but with the updated calibration model reported 
in Sayers et al. (2012) and some other minor modifica- 
tions. 

The key steps involved in our Bolocam data reduction 
and cluster modeling are summarized as follows: We first 
remove sky noise from the time streams by subtracting 
a template of the correlated signal over the field of view 
followed by a high-pass filter (Sayers et al. 2011). This 
process results in a filtered image of the true SZE sig- 
nal (see the left panel in Figure 16), where the filter- 



ing is weakly dependent on the cluster shape due to the 
correlated template removal. To characterize this filter- 
ing, we process a beam-smoothed, initial best-fit cluster 
profile by reverse-mapping it using our pointing infor- 
mation. These data are then processed itcrativcly with 
a new best-fit profile using our full reduction pipeline, 
until the procedure converges. For this analysis we use 
the A10 "universal pressure profile" which adopts a form 
of the Nagai et al. (2007) pressure profile with its slopes 
fixed to the values given in A10, allowing the overall nor- 
malization and scale radius to vary. 

We have derived cluster mass estimates from our SZE 
data alone using the method outlined in Mroczkowski 
(2011). The key innovation of this method is that, in ad- 
dition to assuming HSE, the virial theorem is used, which 
is no stronger an assumption than HSE and can be de- 
rived from HSE and thermodynamics. This method de- 
termines the underlying total mass profile from an SZE- 
determined pressure profile, with the added assumption 
of a constant gas mass fraction / gas . Cluster mass es- 
timates derived with this method have been shown to 
be consistent with X-ray derived results using data from 
the Sunyaev-Zel'dovich Array (SZA, Mroczkowski 2011), 
and SZA followup of blind SZE detections using the At- 
acama Cosmology Telescope (ACT, Reese et al. 2011). 
The SZE-only mass estimates for MACS1206 are given 
in Table 11, which presents Ma and ?*a values at over- 
densities A = 2500 and 500 derived from our Bolocam 
data alone, under the assumptions made. This table also 
contains Bolocam-derived mass estimates at the lensing- 
derived values of r25oo and rsoo- We note that the values 
in Table 11 include an estimate of our systematic errors 
on the SZE-derived masses, which we describe in detail 
below. 

The dominant source of uncertainty in our mass esti- 
mates, as discussed in Mroczkowski (2011), stems from 
the uncertainty in the assumed value of a radially- 
constant f g as(r). Masses derived under this assumption 

— 1/2 

scale as oc / gas ■ We adopt the value / gas = 0.13, 
and marginalize over uncertainties for a range / gas = 
[0.1,0.17], consistent with our X-ray determined gas frac- 
tion measurements at radii near r25oo (see Figure 15). 
An additional source of systematic uncertainty is the ab- 
solute calibration of the Bolocam maps, which is about 
5% and results in a 5% uncertainty in our derived 
masses. Finally, we include a ±1.5% systematic at T25ooi 
and ±5% systematic at rsoo, due to our particular choice 
of parameterization for the pressure profile. These values 
arc roughly consistent with those shown in Mroczkowski 
(2011) for different parameterization of the exponents in 
the pressure profiles. 

By comparison to the lensing mass estimates, we find 
an SZE to lensing mass ratio of 02500 = 1.08±0.29±0.22 
(statistical followed by systematic at 68% confidence) at 
the lensing-derived overdensity radius r25oo of 0.60 Mpc 
(Table 11). Hence, our lensing mass estimate is in agree- 
ment with both the X-ray and SZE mass estimates at 
r 2500- At a lower overdensity of A = 500, we find an 
SZE to lensing mass ratio of 0500 = 1-55 ± 0.30 ± 0.26 
at the lensing derived radius rsoo of 1.3 Mpc, roughly 
consistent with unity within large errors. 

7.4. Dynamical and Physical Conditions of the Cluster 
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MACS 1206 is an X-ray luminous cluster at a redshift 
of z = 0.439, or a cosmic time of t ~ 9 Gyr. The cluster 
appears relatively relaxed in projection in both optical 
and X-ray images, with a pronounced X-ray peak at the 
BCG position (Ebeling et al. 2009; Postman et al. 2012). 
This cluster was classified to be relaxed by Gilmour et al. 
(2009) on the basis of a visual examination of its X-ray 
morphology. Our detailed morphology analysis shows 
no sign of significant recent merging activity around the 
BCG, which is also supported by our strong-lensing anal- 
ysis, finding no significant offset between the DM center 
of mass, BCG, and X-ray peak (Section 4.2). A good 
agreement between the lensing and X-ray mass estimates 
(Section 7.2) indicates that the hot gas is not far from a 
state of hydrostatic equilibrium in the cluster potential 
well. 

However, some evidence of merger activity along the 
line of sight was suggested by the high velocity disper- 
sion of a v rj 1580 km s -1 based on 38 redshift measure- 
ments (Ebeling et al. 2009). Recently, a much larger 
spectroscopic sample of cluster members has been ob- 
tained for this cluster using VLT/VIMOS (P. Rosati et 
al., in preparation). Defining membership is crucial for a 
dynamical analysis since interlopers by projection effects 
can largely bias the derived projected velocity disper- 
sion, especially at large radii where the number density 
of cluster members is low (Wojtak et al. 2007). Using 
a secure sample of > 400 cluster members identified in 
the projected phase space (e.g., Biviano & Salucci 2006; 
Lemze et al. 2009) the velocity dispersion profile is found 
to decrease outward fairly rapidly from ~ 1500 kms" 1 in 
the central region to ~ 800kms _1 at a projected dis- 
tance of R ~ 2 Mpc. Accordingly, the dynamical mass 
estimate is in agreement with the lensing estimate (A. 
Biviano et al., in preparation). This may argue against 
a strong deviation from dynamical relaxation. 

The present Chandra analysis yields a gas temperature 
of 10.8±0.7keV averaged in the radial range [70, 700] kpc. 
Assuming that the galaxies and the gas are confined in 
the same gravitational potential well, this is consistent 
with a line-of-sight velocity dispersion of ~ 1300 km s" 1 , 
which is again in agreement with the observed value. 
This may also suggest that the cluster is not far from 
equilibrium. The cluster appears fairly round in both 
Chandra and XMM images at large distances from the 
cluster center, as demonstrated in Figure 10. X-ray emis- 
sion is concentrated around and peaked on the BCG, 
but shows some elongation within 9 1' at a position 
angle around 120 deg east of north, aligned with the ori- 
entation of the projected mass distribution. The surface 
brightness profile is fairly smooth; but there might be 
some tiny hints of discontinuities (see the ACCEPT cat- 
alog). 50 However, a much deeper observation is required 
to confirm them. 

Finally, morphological analysis of Bolocam data has 
been performed in an identical way to the procedure used 
in Sayers et al. (2011). We find an cllipticity of (10 ± 
7)% with a position angle of 55 ± 27 deg north of west 
from elliptical A10 model fits to our Bolocam SZE data. 
The fits include all data within a 14 x 14 arcmin square, 
corresponding to a fairly large aperture of 9 max 9' > 

50 http:/ /www. pa. msu.edu/astro/MC2/accept/clusters/3277. html 



#vir ~ 7'. Of the approximately 50 clusters observed with 
Bolocam and fit with an elliptical A10 model, MACS1206 
is one of the more circularly symmetric model fits. 

7.5. Comparison with ACDM Predictions 

In Figure 17, we summarize our full lensing con- 
straints on the mass and concentration parame- 
ters of MACS1206, along with recent ACDM pre- 
dictions for relaxed cluster-sized halos based on N- 
body simulations (Duffy et al. 2008; Klypin et al. 2011; 
Bhattacharya et al. 2011; Prada et al. 2011). Our range 
of allowed concentration values (4.6 < c v - u < 7.9 at 
la; see Section 6.6) span the high end and average 
expectations (4 ^ (c v - 1T ) ^ 7) from ACDM simulations 
(Duffy et al. 2008; Zhao et al. 2009; Klypin et al. 2011; 
Bhattacharya et al. 2011; Prada et al. 2011). Average 
concentrations for relaxed clusters are found to be ~ 10% 
higher and have lower scatter than those for the full pop- 
ulation of halos (Duffy et al. 2008; Bhattacharya et al. 
2011). A relatively high concentration of MACS1206 
may also be indicated by the large Einstein radius 
Ei „ « 28"(17") at z s = 2.5 (1.0) (Ebeling et al. 2009; 
Zitrin et al. 2012). 

Care must be taken when comparing these predictions 
for spherically-averaged halo structure with our lensing 
results, which are obtained from an NFW fit to the pro- 
jected lensing measurements assuming a spherical halo. 
In the previous subsection (Section 7.1), we have shown 
that our lensing results are in good agreement with the 
X-ray derived mass profiles (see Figures 14 and 15) in the 
region of overlap ( IMpc), as well as with the Bolo- 
cam SZE mass estimates (Section 7.3), suggesting that 
the level of orientation bias (see Section 7.1.2) is not sig- 
nificant in this cluster. 

Additionally, the effects of baryonic physics can im- 
pact the inner halo profile (at r 0.0 5r v ; r , Duffy et al. 
2010), and thus modify the gravity-only c-M relation, es- 
pecially for less massive halos (M v ; r ^4x 10 14 M Q ft, -1 ; 
see Bhattacharya et al. 2011). Using cosmological hy- 
drodynamical simulations including the back reaction of 
baryons on DM, Duffy et al. (2010) found a < 20% in- 
crease in the halo concentration for cluster-sized halos 
(M vir < 6 x 10 14 M Q h~ l at z = 0). When excluding the 
central 50 kpc /i -1 (« 0.03 r V j r ) region from our primary 
strong-lensing mass model (Zitrin et al. 2012), we find a 
~ 17% decrease in the best-fit concentration parameter 
derived from our full-lensing analysis (Section 6.6), as 
demonstrated in Figure 17. We note, the CLASH clus- 
ters arc massive (5 x 10 14 < M vir /M Q < 3 x 10 15 , sec 
Postman ct al. 2012), and hence expected to be less af- 
fected by baryonic effects. 

For this cluster, the lensing-derived total mass distri- 
bution is consistent with the NFW form (a = 73d(?" — > 
0) = 0.961q'4q), as found for several relaxed clusters: 
A611 (Newman ct al. 2009); A383 (Zitrin et al. 2011c); 
A1703 (a w 0.9 Richard et al. 2009; Oguri et al. 2009); a 
stacked full-lensing analysis of A1689, A1703, A370, and 
C10024+17 (a = 0.89^;^, Umetsu et al. 2011a). Mul- 
tiwavelength observations can be used to measure gas 
and stellar density profiles for subtraction from lensing- 
derived total mass profiles to yield DM-only mass pro- 
files (Lemze et al. 2008; Newman et al. 2009), allowing 
for a more direct comparison with CDM predictions from 
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gravity-only simulations, 
forthcoming paper. 



We defer this analysis to a 
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Figure 17. Constraints on the halo mass and concentration 
parameters (M v j r , c v i r ) for the X-ray selected CLASH cluster 
MACS1206 (2 = 0.439) derived from spherical NFW fits to com- 
bined weak and strong lensing observations, compared to ACDM 
predictions for relaxed populations of simulated cluster-sized halos 
at 2 = 0.44 (except 2 = 0.5 for Klypin et al. 2011). Our results 
are shown with and without the central 50kpc?i~ 1 cut (Section 
6.6) applied to the Zitrin et al. (2012) based strong-lensing model. 
The A^-body predictions by Duffy et al. (2008) and Prada et al. 
(2011) are shown in light blue, including lcr lognormal scatter 



(0.11 in log! 



29%) indicated by the respective hatched 



areas. Portions of these lines are dashed to indicate extrapolations 
to higher masses. Average results from two additional simulations 
(Klypin et al. 2011; Bhattacharya et al. 2011) are shown in blue 
for clarity. Duffy et al. (2008) and Bhattacharya et al. (2011) de- 
rived results for dynamically relaxed cluster subsamples, yielding 
concentrations ~ 10% higher than for the full samples. This 10% 
factor has been applied to the results from the other simulations. 



8. SUMMARY AND CONCLUSIONS 

In this paper, we have presented a comprehensive lens- 
ing analysis, combining independent measurements of 
wcak-lcnsing distortion, magnification, and strong lens- 
ing of the massive X-ray selected cluster MACS1206 at 
z = 0.439. This is based on wide-field Subaru BVR c I c z' 
imaging, combined with detailed strong-lensing informa- 
tion obtained from deep CLASH HST 16-band imaging 
and VLT/VIMOS spectroscopy (Zitrin et al. 2012). 

The deep Subaru multi-band photometry is used to 
separate background, foreground and cluster galaxy pop- 
ulations using the selection techniques established in 
our earlier work (Medezinski et al. 2010; Umctsu et al. 
2010), allowing us to obtain a reliable weak lensing sig- 
nal free from significant contamination of unlcnscd clus- 
ter and foreground galaxies. By combining complemen- 
tary distortion and magnification measurements, we con- 
structed a model-free mass distribution out to well be- 
yond the virial radius (r v - lr w 1.6 Mpc h^ 1 ). In addi- 
tion to breaking the mass-sheet degeneracy inherent in 
shape distortion measurements, the magnification mea- 
surements also increase the overall significance by ~ 23% 
(Section 5.3). 

We have also obtained an improved inner mass distri- 
bution from a reanalysis of the Zitrin ct al. (2012) data 
using our new MCMC implementation of the Zitrin et al. 



(2009) method. We introduced a technique to self- 
calibrate the bin-to-bin covariance matrix of the inner 
mass profile (Section 4.3), accounting for possible sys- 
tematic errors inherent in the analysis. This is a cru- 
cial step for a joint analysis to combine constraints in 
different regimes of signal strength. The inner radial 
boundary for the mass profile is chosen to be sufficiently 
large to avoid smoothing from cluster miscentering effects 
(Johnston et al. 2007). The derived inner mass profile is 
shown to be consistent with our semi-independent re- 
sults from a wide variety of four strong-lensing analyses 
(Lenstool, Pixelens, LensPerfect, and SaWLens; 
see Section 4.4), and overlap well with the Subaru-based 
outer mass profile, ensuring consistency in both the weak 
and strong regime. 

The Subaru data reveal the presence of an elongated 
large scale structure around the cluster, in both the dis- 
tribution of galaxies and from the mass distribution, 
with the major axis running NW-SE, aligned well with 
the cluster and BCG shapes, showing elongation with a 
^2:1 axis ratio in the plane of the sky (Section 6.2). 
The azimuthally-averagcd projected mass profile from 
our full-lensing analysis exhibits a shallow profile slope 
dlnE/dlni? ~ -1 at cluster outskirts (R > IMpc/i" 1 ), 
whereas the mass distribution excluding the NW-SE ex- 
cess regions steepens further out, well described by the 
standard NFW form (Section 6.4). Assuming a spherical 
halo, we have obtained a virial mass M v ; r /1O 15 M0 h^ 1 = 
1.07lg '^g(stat.) ± 0.10(syst.) and a halo concentration 
Cyir = 6.9lg'g(stat.) ± 1.2(syst.), which is somewhat high 

but falls in the range 4 ^ (c) ,$ 7 of average c(M, z) pre- 
dictions for relaxed clusters from recent ACDM simula- 
tions. When the innermost 50kpc/i _1 is excluded from 
the fit, we find a slightly lower concentration c v j r = 
5.Tti'*(stat.) ± 1.2(syst.), a decrease of approximately 
17% (Section 6.6). 

We have shown that our full-lensing mass profile is in 
agreement with Chandra+XMM X-ray data in the region 
of overlap (Figure 14). The hydrostatic X-ray to lensing 
mass ratio, oa = Mx(< r^)/Mi ons (< r&), is consistent 

with unity especially at r ~ r25oo with 02500 = 0.95^ ; 2 5- 
Our full lensing results, when combined with Chandra gas 
mass measurements, yield a gas mass fraction estimate 
free from the HSE assumption. We find a cumulative gas 

mass fraction of / gas (< t) = 13.7l" 3 ;o% at r ~ l-7?"2500) 
a typical value observed for high mass clusters. Overall 
good agreement is also obtained with SZE-only cluster 
mass estimates based on Bolocam observations (Section 
7.3). 

The CLASH survey is producing substantial improve- 
ments in both the quality and quantity of direct empiri- 
cal constraints on cluster-sized DM halos (Postman et al. 
2012; Zitrin et al. 2011c; Coe et al. 2012; Zheng et al. 
2012), for an X-ray selected sample of relaxed clusters, 
selected free of lensing selection bias, as well as for a 
lcnsing-selected sample of high-magnification clusters. 
The CLASH imaging, in combination with Subaru weak- 
lcnsing data, allows us to make precise measurements of 
the mass distributions of individual clusters over the full 
range of cluster radius, and to help understand the possi- 
ble evolutionary and tidal effects of connecting filaments 
and local clusters on the mass distribution of the cen- 
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tral cluster, for a detailed comparison with the standard 
ACDM cosmology and a wider examination of alterna- 
tive scenarios. With the full sample of CLASH clusters, 
we will be able to establish the representative mass pro- 
file of massive clusters in gravitational equilibrium, and 
robustly test models of structure formation. 
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Table 1 

Properties of the galaxy cluster MACS1206 



Parameter 



ID 

Optical center position (J2000.0) 

R.A 

Decl 

X-ray peak position (J2000.0) 

R.A 

Decl 

Redshift 

X-ray temperature (keV) 

Einstein radius (") 



Value 



MACS J1206.2-0847 

12:06:12.15 
-08:48:03.4 

12:06:12.28 
-08:48:02.4 

0.4385 
10.9 ±0.6 
28 ± 3 (17 ± 2) at z 3 = 2.54 (1.03) 



References. — [1] Ebeling et al. (2009); [2] Zitrin et al. (2012); [3] 
Postman et al. (2012); 

Note. — The cluster MACS J1206. 2-0847 (z = 0.4385) was discov- 
ered in the Massive Cluster Survey (MACS) as described by Reference 
[1] . The optical cluster center is defined as the center of the BCG from 
Reference [2] . Units of right ascension are hours, minutes, and seconds, 
and units of declination are degrees, arcminutes, and arcseconds. The X- 
ray properties are taken from Reference [3]. See also Reference [1]. The 
BCG is located within Ri 2" (a projected separation of 9kpc/i _1 ) of 
the X-ray emission peak. The Einstein radii are constrained by detailed 
strong lens modeling by Reference [2], 



Table 2 

Subaru Suprime-Cam data. 



Filter Exposure time a Seeing* 3 mii m c 

(ks) (arcsec) (AB mag) 



B 2.4 1.01 26.5 

V 2.2 0.95 26.5 

R c 2.9 0.78 26.2 

I c 3.6 0.71 26.0 

z' 1.6 0.58 25.0 



a Total exposure time. 

b Seeing FWHM in the full stack of images. 

c Limiting magnitude for a 3<r detection within a 

2" aperture. 



Table 3 

Galaxy color selection. 



Sample 


Magnitude limits' 1 


N 








(AB mag) 




(arcmin 2 ) 




Red 


21.5 < z' < 24.6 


13252 


9.9 


1.16 


Green 


z' < 24.6 


1638 


3.4 


0.44 


Blue 


22.0 < z' < 24.6 


4570 


4.3 


1.95 



a Magnitude limits for the galaxy sample. 
b Mean surface number density of source background galaxies. 
c Mean photometric redshift of the sample obtained with the 
BPZ code. 



Zitrin, A., Broadhurst, T., Barkana, R., Rephaeli, Y., & Bem'tez, 

N. 2011a, MNRAS, 410, 1939 
Zitrin, A., Broadhurst, T., Coe, D., et al. 2011b, MNRAS, 237 
Zitrin, A., Broadhurst, T., Umetsu, K., et al. 2009, MNRAS, 396, 

1985 



— . 2010, MNRAS, 408, 1916 

Zitrin, A., Broadhurst, T., Coe, D., et al. 2011c, ApJ, 742, 117 
Zitrin, A., Rosati, P., Nonino, M., et al. 2012, arXiv, 1107.2649 
Zwicky, F. 1959, Handbuch der Physik, 53, 390 
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Table 4 

Galaxy samples for weak-lensing shape measurements. 



Sample 


JV 


n 9 a 




2 S 


lC ff c 


{Di s /D s ) d 






(arcmin 2 ) 




M1206 


COSMOS 


M1206 COSMOS 


Red 


8969 


9.2 


0.42 


1.05 


1.05 


0.51 0.51 


Blue 


4154 


4.3 


0.48 


1.55 


1.58 


0.62 0.63 


Blue+red 


13123 


13.4 


0.44 


1.15 


1.12 


0.54 0.53 



Note. — 

a Mean surface number density of source background galaxies. 

b Mean rms error for the shear estimate per galaxy, a g = (cr 2 ) 1//2 . 

c Effective source redshift corresponding to the mean depth (/3) of the sample. 

d Distance ratio averaged over the redshift distribution of the sample, 



Table 5 

Best-fit NFW model parameters for MACS1206. 



M vir a 
(1O 15 M h- 1 ) 


Cvir 


X 2 /dof 


0Ein b 
(") 


N c 
WL 


SL 




Method 


qq+0.32 
u - t,a -0 .26 

1 ic+0.34 
1 - lo -0.28 
i ic+0.25 
i ' ±o -O.20 


c y+3.6 
°-'-2.1 
4 0+ 2 ' 1 

'•°-1.8 


3.3/6 
4.5/7 
10.6/6 


21 
11 

32 


8 
9 
8 







1) WL tangential distortion (§6.1.1)) 

2) WL tangential distortion + magnification (§6.1.1)) 

3) WL(#2) + LSS correction d (§6.4)) 


0-88-0.21 
0.97±g;S 




3.9/8 
6.9/9 


28 
26 


8 
9 


2 
2 


4) WL(#1) - 

5) WL(#2) - 


h Einstein radius (§6.1.2)) 
h Einstein radius (§6.1.2)) 


1 14+ U -^ 

i - ii -0.18 

i n7+ ' 20 


6 - 6 -0.9 

6 9+ 10 
D - M -0.9 


24.2/31 
18.0/31 


28 
28 


8 
8 


25 
25 


6) WL(#2) - 

7) WL(#3) - 


h SL f (§6.5)) 

h SL = our primary NFW result (§6.5)) 



Note. — All our methods take into account the cosmic covariance from distant, uncorrelated large scale structure 
(LSS) projected along the line of sight. For weak lensing, the source redshift uncertainty (^ S)G ff — 1.15 ± 0.1) of our 
background sample has been marginalized over. 



a The virial overdensity is A v i r w 132 times the critical density of the universe at z — 0.439 in the adopted cosmology 

(Q m = o.3, n A = 0.7). 

Effective Einstein radius for a source at z s — 2.5 predicted by the model. The observed value is 28" ± 3". 
c Respective numbers of weak and strong lensing constraints. 

^ Excluding the elongated LSS around the cluster extending along the NW-SE direction (see Figure 8). 

c Combining with double Einstein-radius constraints of #Ein — 17 ±2 at z s — 1.03 and #Ein — 28 ±3 at z s — 2.54. 

Additionally, an rms displacement of 2 'is assumed for each system due to uncorrelated LSS projected along the line of 

sight, and is combined in quadrature with the respective measurement error to estimate a total uncertainty. 

* Combining with the inner strong-lensing based mass profile derived from an MCMC implementation of Zitrin ct al. 

(2012) (Sections 4.1—4.3). The outer fitting radius is limited to less than 12' for direct comparison with Method #7 based 

on the LSS-corrected weak-lensing profile (Section 6.4). 



Table 6 

Ellipticity and position angle measurements. 



Method 


a a 
"max 


Ellipticity b 


PA C 




(') 




(deg) 


BCG 


10" 


n 5 o + 0.03 


l5.0tH 


Chandra X-ray 


1.5' 


sn+ 03 


21.9+i- 7 


Galaxy # density 


8' 


5^+ 04 


15 7+ 13 
10 -'-5.9 


Galaxy light 


8' 


41 + ' 06 
u '^ 1 -0.06 


19.0+J-9 


WL mass map 


8' 


0.371^3 


19.4±f 7 5 7 


WL 2D shear fit 


8' 




28.6l 7 ;g 



a Circular aperture radius. 

^ Ellipticity modulus defined such that, for an ellipse with 

major and minor axes a and b, it reduces to 1 — b/a. 

c Position angle of the major axis measured north of west. 
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Table 7 

Best-fit gencralizcd-NFW model parameters for MACS1206. 



Method 3 


M vir 
(IO^Mq/i" 1 ) 


c_ 


2 b 




X 2 /dof 


#Ein d 
(") 


#6 
#7 


i i 7 +0.29 
1 - xl -0.22 
1 06+ 023 


6.3" 
7.0 


(-1.5 
-1.5 
(-1.5 
-1.4 


1 09+ ' 28 
±.uw_ Q 42 

96+ 031 
u.»o_ 4g 


24.1/30 
18.0/30 


28 
28 



Note. — See for details Section 6.5. 
a Fitting method in Table 5. 

k Effective concentration parameter for gNFW, c„2 = — 
Cv, r /(2 - a). 

c Central cusp slope of gNFW. 

^ Effective Einstein radius for a source at z s — 2.5 predicted by the model. 
The observed value is 28" ± 3". 



Table 8 

Impact of the choice of strong lensing models in the full lcnsing analysis. 



A'/ V ir 

(IO^A/q/i- 1 ) 


( 'vir 


X 2 /dof 


(") 




Method b 


1 07+ ' 20 
i 17 +o.25 
± - ± 1 -0.20 
, 07+0.26 

J "°'— 0.22 

1 26+ ' 20 
i - ZD -0.17 


5.7tH 

c- O+0.9 

& - 8 -0.8 
6 +0 !! 


18.0/31 
16.0/26 
15.4/20 
11.6/31 


28 
25 
29 
28 


WL(#3) - 
WL(#3) - 
WL(#3) - 
WL(#3) - 


h Zitrin+12 c 

f- Zitrin+12 + 50kpc/h cut d 

- PlXELENS 

- Lenstool 



Note. — See for details Section 6.6. 
a Effective Einstein radius for a source at z s — 2.5 predicted by the model. 

Combination of strong and weak lensing mass models used for the fitting. For all cases, 
Method #3 of Table 5 is used for weak lensing. 
c This corresponds to our best model (Model #7) of Table 5. 

Now applying a central 12" (~ 50kpcfr -1 ) cut to the Zitrin et al. (2012) based strong 
lcnsing model. 



Table 9 

Best-fit SPLS model parameters for MACS1206. 



Method 3 




(kpeft- 1 ) 




X 2 /dof 


A/vir 


K 


(kpch- 1 ) 


V 


X 2 /dof 


(10 15 Af Q h- 1 ) 


#6 
#7 


r 7 +0.70 
-0.54 
fi2 +0.71 


23.ll 44 
22.4± 4 ;° 


1 
1 


26.5/31 
23.9/31 


1.78 ±0.56 
1.74 ±0.55 


1 7+O.6O 
-0.48 
07+O.57 
J - u ' -0.47 


33.0ti° 4 9 
36.8t^ 2 4 8 


o.84±°:ii 
o.r7±8:« 


24.0/30 
19.5/30 


1.41 ± 0.34 
1.26 ±0.37 



Note. — The convergence profile of the softened power-law sphere (SPLS) model, p(r) — po(l + r2 /' r c)' T? 3 " 2 5 is given by — 
kq (l + 8 2 /^c)^ where #c — fc/Di and kq — B (.j, 1 — ^) por c /T, CT it with B being the standard Euler beta function. Here X cl -it is 

evaluated for a source at a reference redshift of z s — 2.54. For details, see Section 6.7. 
a Fitting method in Table 5. 



Table 10 

Comparison with X-ray cluster mass estimates. 



Data 


A-/2500 


C2500 


^2500 


a a 
y Ein 




(10 14 M Q ) 




(Mpc) 


(") 


Chandra 


4± 1 


1.8 ± 1.5 


0.6 ± 0.1 


23 


Chandra+XMM 


4.5 ±0.3 


0.9 ±0.3 


0.58 ±0.02 


20 


WL+SL b 


4.9 ± 1.3 


1.8 ±0.3 


0.60 ±0.06 
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Note. — See for details Section 7.2. All quantities here are 
given in physical units assuming the concordance ACDM cosmology 
(h = 0.7, flm = 0.3, n A = 0.7). 



a Effective Einstein radius for a source at z s — 2.5 predicted by the 
model. 

Model #7 of Table 5 based on the full weak and strong lensing 
constraints. 



Umetsu et al. 
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Table 11 

Bolocam SZE-dcrived cluster mass estimates. 



Overdensity a 
A 


Bolocam-derived r A b 

r A M(<r A ) 
(Mpc) (10 14 M o ) 


Lensing-derived r A c 
r A M(< r A ) 
(Mpc) (10 14 M Q ) 


2500 
500 


n K „+0.01+0.06 
UitK >-0.02-0.05 
i fi7 +0.09+0.12 
-0.08-0.12 


r O+0.4+1.7 
°'°-0.4-1.4 
91 9+3.7+5.1 
zl -3.0-4.3 


0.60 
1.31 


c .,+0.2+0.8 
3-J_0.2-0.7 
i c 7 +1.2+2.3 
10 -'-l. 1-2.1 



Note. — For each value the first error estimate represents our measure- 
ment uncertainty and the second error estimate represents our uncertainty 
due to systcmatics in our fitting method, flux calibration, and choice of pa- 
rameterization. Sec for details Section 7.3. All quantities here arc given in 
physical units assuming the concordance ACDM cosmology (h — 0.7, n m — 

0.3, n A = 0.7). 

a Mean interior ovordensity with respect to the critical density of the universe 
at the cluster redshift z — 0.439. 

Bolocam cluster mass estimates at the Bolocam-SZE derived values of over- 
density radius 

c Bolocam cluster mass estimates at the lensing-derived values of overdensity 
radius ta ■ 



APPENDIX 

A. COMBINING LENS DISTORTION AND MAGNIFICATION 
A.l. One- Dimensional Method 

We first summarize the Bayesian method of Umetsu et al. (2011b) for a direct reconstruction of the cluster mass 
profile from combined radial distortion and magnification profiles. 

In the Bayesian framework, we sample from the posterior PDF of the underlying signal s given the data d, P(s\d). 
Expectation values of any statistic of the signal s shall converge to the expectation values of the a posteriori marginal- 
ized PDF, P(s\d). For a mass profile analysis, s is a vector containing the discrete convergence profile, K, = 
with i = 1,2, .., N in the subcritical regime > #Ein)j an d the average convergence within the inner radial boundary 
0min of the weak lensing data, K nnn = «(< # m in), so that s = {7S m j n , Ki}fL 1: being specified by (N + 1) parameters. 

Baycs' theorem states that 

P(s\d) oc P(s)P(d\s), (Al) 

where C(s) = P(d\s) is the likelihood of the data given the model (s), and P(s) is the prior probability distribution 
for the model parameters. The C(s) function for combined weak lensing observations is given as a product of the 
two separate likelihoods, C = C g+ C lll where C g+ and £ M are the likelihood functions for tangential distortion and 
magnification bias, respectively. The log-likelihood functions for the weak-lensing observations {g+,i}f = i and {n^^i}f =l 
are given respectively (ignoring constant terms) as 

W^-ln^^EN^, (A2) 
z i=i 

Us) - - ln£, = \ £ (A3) 

Z i=l °*M 

where {g+,i}f = \ and {n fJli i}fL 1 are the theoretical predictions for the corresponding observations. The total likelihood 
l\Y)(s) = — ln£ of the combined observations is obtained as 

hr>=l g+ +l^ (A4) 

Here we consider a simple flat prior with a lower bound of s = 0. Additionally, we account for the uncertainty in the 
calibration parameters, c = (no, s, ui), namely the normalization and slope parameters (no, s) of the background counts 
and the relative lensing depth lo = (ft (red) )/(/3 (back)) between the background samples used for the magnification 
and distortion measurements. 

We use the Markov chain Monte Carlo (MCMC) technique with Metropolis-Hastings sampling to constrain our mass 
model s. The covariance matrix C of s is obtained from MCMC samples. 

A. 2. Two- Dimensional Method 

Here we extend the ID method of Umetsu et al. (2011b) to a 2D mass distribution n(8), by combining 2D distortion 
data with the azimuthally-averaged magnification information. For this analysis, the signal s is a vector of parameters 
containing discrete mass elements on a 2D Cartesian grid of independent cells: s = {K m }^"". The 7(6) field can be 
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written as a linear combination of the parameters s (Equation (3)). Then, the distortion g(0) and magnification fj,(8) 
fields can be uniquely specified in the subcritical regime (Section 2). 

In analogy to Equation (13), we calculate the weighted average g a , m = 9a(9m) {ol = 1,2) of individual distortion 
estimates, and its covariance matrix, 

Cov[g a , m ,gp,n\ = {C a ) af}mn = ^<r 2 g (9 m )S m „S a p, (A5) 

where a g (9 m ) is the standard error of the weighted mean distortion, g(6 m ). Accordingly, the 2D shear log-likelihood 
function l g (s) = — ln£ g is written as 

l N cli u 2 

l g( S ) = 2 Yl 12 ^ 9a ' m ~ 9^m( s )} (Wj) Q/i) „ [9ff,n - 9p,n( 8 )]i ( A6 ) 
m,n— 1 a, j3= 1 

where g a ,m(s) is the theoretical prediction for #Q. m , and (yV g ) a /3.mn is the shear weight matrix, 

(W fl ) a/s>mn = M m M n (C fl - 1 ) a/}fTOn , (A7) 

with M m being a mask weight, defined such that M m = if the mth cell is masked out and M m = 1 otherwise. 
In practice, we exclude from our analysis innermost cells which lie in the cluster central region, where the surface 
mass density can be close to or greater than the critical value (i.e., k ^ 1). Furthermore, this is crucial to minimize 
contamination by unlensed cluster member galaxies (see Section 3.4). 
Now we combine 2D distortion data with magnification information to obtain the total log- likelihood ?2d( s ) as 

ho =l g + ^, (A8) 

where L, given by Equation (A3), imposes a set of azimuthally- integrated constraints on the underlying k field. 
Since the degree of magnification is locally related to k, this will essentially provide the (otherwise unconstrained) 
normalization of n(6) over a set of concentric rings where count measurements n^^ are available. Note, no assumption 
is made of azimuthal symmetry or isotropy of the cluster mass distribution. 

This 2D inversion problem involves estimation of a large number of parameters s; typically, N ce \\ ^ 1000 when 
distortion data are binned into subarcminute pixels. We use in our implementation the conjugate-gradient method 
(Press et al. 1992) to find the best solution. We include Gaussian priors on the calibration nuisance parameters 
c = (s,tiq,uj), given by means of quadratic penalty terms with mean values and variances directly estimated from 
data. The log posterior PDF, F = — h\P(s\d), is expressed as a linear sum of ^d(s) and the prior terms on c. 
The best-fit parameters are determined with a maximum-likelihood estimation, by minimizing the function F with 
respect to p = (s, c), a vector containing the mass and calibration parameters. Here we employ an analytic expression 
for the gradient function VF(p) obtained in the nonlinear subcritical regime. To quantify the errors on the mass 
reconstruction, we evaluate the Fisher matrix at the maximum likelihood estimate p = p as 



, d 2 F(p) 

J~ r. 



dp m dpn 



(A9) 

p=p 



where the angular brackets represent an ensemble average, and the indices (m, n) run over all model parameters. We 
estimate the covariance matrix C of s as 

Cmn = i?- 1 )^ ■ (A10) 



